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EDITORIAL 

The  direct  time-domain  modelling  of  electromagnetic  fields  and  high-frequency  circuits  meets  with 
growing  interest.  Modem  powerful  computers  make  feasible  the  applications  of  time-domain 
methods  in  the  modelling  of  electromagnetic  fields  and  networks.  The  advantages  of  time-domain 
methods  are  their  high  flexibility,  their  potential  to  include  non-linear  effects  and  time-dependent 
parameters,  and  their  transparency  with  respect  to  concepts  and  algorithms.  Time-domain  analysis 
elucidates  the  physical  principles  underlying  the  phenomena  and  supports  a  creative  design  of 
circuits  and  systems.  For  these  reasons,  time-domain  methods  are  of  high  interest  for  the  develop¬ 
ment  of  CAD  tools  for  the  modelling  of  microwave  and  millimetre-wave  integrated  circuits,  and 
broad-band  microwave  devices,  antennas,  circuits  and  systems.  The  combination  of  field  concepts 
and  network  concepts  allows  the  segmentation  of  complex  structures  and  to  apply  full  wave 
analysis  to  the  segments. 

This  special  issue  is  the  second  of  three  parts  comprising  contributions  to  the  workshop  on  the 
German  IEEE  MTT/AP  Joint  Chapter  and  the  German  IEEE  CAS  Chapter  on  Discr^  Time 
Domain  Modelling  of  Electromagnetic  Fields  and  Networks  on  24  and  25  October  1991  at  the 
Technische  Universitat  Miinchen.  The  first  part  was  published  in  vol.  5.  no.  3  of  the  Journal.  The 
purpose  of  this  workshop,  organized  by  Peter  Russer  and  Josef  Nossek  under  the  sponsorship  of 
the  European  Research  Office  (ERO)  of  the  US  Army,  was  to  bring  together  researchers  dealing 
with  time -domain  simulation  and  transient  phenomena  in  fields  and  networks. 


Peter  Russer 


Hi. 
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EFFICIENT  ANALYTICAL-NUMERICAL  MODELLING  OF 
ULTRA-WIDEBAND  PULSED  PLANE  WAVE  SCATTERING 
FROM  A  LARGE  STRIP  GRATING 

LAWRENCE  CARIN  AND  LEOPOLD  B.  FELSEN 

Weber  Research  Institute/ Electrical  Engineering  Department,  Polytechnic  University.  Farmingdate,  NY  11735,  U.S.A. 


SUMMARY 

Ultra-wideband  (UWB)  pulsed  plane  wave  scattering  from  a  large  but  finite  stnp  grating  in  free  space  is 
analysed  in  the  frequency  domain  via  decomposition  into  plane  wave  spectra,  implemented  numerically  by 
the  method  of  nionients,  and  then  inverted  into  the  time  domain  (TD).  To  make  this  procedure  practical 
under  UWB  conditions,  closed  form  expressions  are  derived  for  interaction  integrals  involving  widely 
separated  expansion  and  testing  functions.  These  closed  forms  are  based  on  a  judicious  choice  of  the  basis 
functions,  and  on  asymptotic  methods  for  reducing  the  integrals.  Although  large  separation  distances  are 
assumed,  the  expressions  have  been  found  to  be  accurate  for  separations  as  small  as  01  wavelengths.  The 
TD  self  terms  can  also  be  evaluated  efficiently.  To  test  the  frequency  domain  aigonthm.  compansons  are 
made  with  available  data  in  the  literature  for  surface  currents  and  far-field  scattering  from  multiple  strips. 
New  short  pulse  TD  results  are  shown  as  well- 


1.  INTRODUCTION 

Plane  wave  scattering  from  a  collection  of  periodically  arranged  elements  continues  to  be  a  topic 
of  interest.  Periodic  arrays  of  patches  or  slots  have  been  used  for  microwave  and  millimetre-wave 
frequency  selective  surfaces.'  Strip  gratings  have  found  use  in  optical  spectrometers^  and  as 
dispersive  elements  in  pulse  compression  systems.'  Although  truncated  in  space,  the  arrays  are 
usually  electrically  large  and  are  therefore  often  treated  by  analysing  an  ideal  infinitely  periodic 
array.  In  such  studies,  the  problem  reduces  to  the  much  simpler  investigation  of  scattering 
from  a  single  unit  cell.  Recently,  however  attention  has  been  given  to  the  effects  of  array 
truncation.'-^' 

Nearly  all  investigations  of  scattering  from  arrays  of  elements  have  been  performed  in  the 
frequency  domain. With  current  interest  in  impulse  or  UWB  radar."  the  time-dependent 
scattering  of  short  pulses  from  such  configurations  is  gaming  in  importance.  Moreover,  the 
availability  of  picosecond  and  femtosecond  lasers  makes  these  studies  relevant  also  the  interac¬ 
tion  of  infrared  or  optical  pulses  with  gratings.  It  is  the  purpose  of  this  paper  to  develop  an 
efficient  technique  for  the  analysis  and  numerical  calculation  of  UWB  pulse  scattering  from  a 
large  but  finite  collection  of  elements.  The  basic  phenomena  associated  with  such  scattering  can 
be  modelled  by  the  strip  array  prototype  adopted  here. 

For  UWB  radars,  the  commonly  accepted  definition  of  a  UWB  pulse  is  one  having  a  bandwidth 
of  25  per  cent  or  more  with  respect  to  the  centre  frequency.®  For  the  present  study,  an  alternative 
definition  is  more  appropriate;  the  UWB  pulse  must  contain  sufficient  energy  at  wavelengths 
ranging  from  Kn  <  D  to  Ko>  D,  where  D  is  the  characteristic  size  of  the  scatterer;  this  range  of 
wavelength  accommodates  at  the  extremes  high  resolution  of  local  features  as  well  as  collective 
wave  phenomena  associated  with  global  features.  To  develop  techniques  for  the  general  analysis 
of  UWB  scattering  from  a  large  but  finite  collection  of  ele.nents,  an  array  of  planar  strips  in  free 
space  has  been  selected  as  a  prototype  problem.  For  this  case,  which  is  of  interest  in  its  own  right, 
the  characteristic  size  D  for  the  UWB  pulse  is  the  strip  width. 

To  analyse  UWB  scattering  efficiently,  special  considerations  must  be  addressed.  If  the  problem 
is  first  analysed  in  the  frequency  domain  and  then  converted  to  the  time  domain  via  the  Fourier 
transform,  thousands  of  frequency  points  are  often  required  to  get  accurate  time-domain  results. 
If  one  were  to  apply  previously  developed  frequency-domain  techniques' directly  to  such  a 
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problem,  the  CPU  time  required  would  be  so  excessive  as  lo  make  the  analysis  impracticable.  To 
avoid  this  difficulty,  the  present  study  utilizes  a  hybrid  numerical/analytical  technique.  This 
involves  application  of  a  spectral-domain  formulation,  with  a  moment  method  solution.  Closed 
form  asymptotic  expressions  are  developed  for  reaction  integrals  that  contain  expansion  and  testing 
functions  separated  by  O  lXo  or  more.  This  method  leads  to  a  highly  efficient  and  accurate 
procedure. 

The  paper  is  organized  as  follows.  Section  2  deals  briefly  with  the  spectral  domain  formulation 
of  time-harmonic  plane  wave  scattering  from  a  collection  of  strips  in  free  space.  Sections  3  and  4 
are  concerned  with  the  techniques  proposed  to  make  such  a  formulation  practicable  for  UWB 
pulsed  scattering  applications.  In  particular,  the  basts  functions  and  integration  techniques  are 
discussed  in  detail.  Numerical  results  are  presented  in  Section  5.  Comparisons  are  made  with 
available  frequency-domain  data  in  the  literature,  followed  by  new  time-domain  results.  Con¬ 
clusions  that  can  be  drawn  from  this  work  are  summarized  in  section  VI. 


2.  FORMULATION  AND  FREQUENCY-DOMAIN  SOLUTION  STRATEGY 

This  section  deals  with  time-harmonic  plane  wave  scattering  from  a  finite  array  of  perfectly 
conducting  infinitesimally  thin  strips  in  free  space.  Referring  to  Figure  1.  the  surfaces  of  the 
various  strips  are  assumed  to  be  perpendicular  to  y.  and  the  fields  in  this  two-dimensional  problem 
are  assumed  to  be  independent  of  z.  Unlike  previous  studies  that  have  performed  the  analysis  by- 
using  the  two-dimensional  free  space  (space  domain)  Green  s  function. the  problem  is  formu¬ 
lated  here  in  the  spectral  domain  (with  respect  to  .r).  This  is  done  for  two  reasons;  (i)  as  shown 
in  section  3.  one  obtains  thereby  a  convenient  and  efficient  asymptotic  representation,  and  (ii) 
this  method  is  readily  extended  to  more  complicated  configurations  involving  layered  dielectrics. 
Because  spectral  domain  formulations  have  been  used  for  several  related  problems.''’'"  this  section 
contains  only  a  brief  summary  of  those  issues  which  are  of  importance  for  the  present  investigation. 

Assuming  a  plane  wave  incident  obliquely  on  the  strips  in  Figure  1 .  and  applying  the  boundary 
condition  for  the  electric  field  on  the  perfect  conductors,  one  arrives  at  the  expression; 

yx(E'-fE’)  =  0  (1) 

A  bold-face  symbol  denotes  a  vector  quantity  and  a  tilde,  later  on.  identifies  quantities  in  the 
spectral  domain.  The  boundary  condition  in  (1)  is  applied  on  the  surface  of  each  strip.  Here. 
E'(x,y)  is  the  incident  vector  electric  field  in  the  absence  of  the  strips,  while  E'(x.y)  is  the  scattered 
vector  electric  field  produced  by  the  electric  surface  currents  J(x'.>')  induced  on  the  strips.  The 
scattered  field  can  be  expressed  as  (an  e'""  time-dependence  is  assumed  and  suppressed  henceforth) 

E’(jr,y)  =  |^G(k„y;y')j(k„y')e-'*.''-‘'dA:,  (2) 


Figure  ).  Example  of  a  multilayer  strip  grating  in  free  space 
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where  G{kj„y;y')  and  are  the  dyadic  Green's  function  and  surface  current,  respectively, 

in  the  it,  spectral  wavenumber  domain.  The  contour  of  integration  C  is  assumed  to  run  initially 
along  the  real  axis.  Because  the  problem  is  two-dimensional,  only  a  single  component  of  surface 
current  (longitudinal  or  transverse)  is  induced  for  a  given  polarization  (TE  or  TM.'*  respectively). 
Therefore,  only  a  single  component  of  the  dyadic  Green  s  function  is  required  for  a  given  incident 
polarization.  The  required  component  of  the  spectral  domain  Green's  function  is 


GHk„y,y') 


Iv-v  I 


(3) 


with  ko  =  Fot  TM  and  TE  waves,  respectively,  it  is  well  known  that  the  wave  impiedances 

are  given  by 


Z”  = 


\ki-kr, 

jaiCfl 


\.kl-kr, 


(4) 


The  boundary  condition  in  (1)  is  enforced  numerically  by  first  expanding  the  unknown  surface 
currents  Jfx'.y')  in  a  known  set  of  basis  functions  ft(jr'.y  )  with  unknown  cocfticients  a*: 

=  E (5) 
By  applying  a  Galerkin  testing  procedure*'  to  (1).  one  obtains 

-fE-(jc.y)-f*(x,y)(k=  ^  f  fA(M')-G(^,.y;>')-^(*,.>  )e  '‘.'v-.- u' d*,  (6) 

Js  -  1  jr 

for  m  =  1,  2,  ....  N^.  Here,  x,  represents  the  location  (in  x)  of  the  centre  of  the  basis  function 
f„  and  superscript  *  denotes  the  complex  conjugate.  The  integral  on  the  left  side  of  (6)  extends 
over  the  surface  .S'  of  a  particular  strip.  By  expanding  the  currents  and  applying  the  testing 
procedure  on  each  strip,  an  Nf,xNt,  matrix  equation  is  produced,  where  .V^  is  the  total  number 
of  basis  functions.  From  this  equation  one  can  determine  the  basis  function  coefficients,  and  from 
(5)  and  (2),  respectively,  the  currents  and  scattered  fields. 

For  a  plane  wave  incident  at  angle  0^  (see  Figure  1),  E‘(x..v)  can  be  expressed  as 


E'(x,y)  =  c  e'*o*'"  (7) 

where  c  is  a  vector  constant.  For  TE  incidence,  c  is  in  the  z-direction  while  for  TM  incidence,  c 
lies  in  the  x-y  plane  but  depends  on  the  angle  of  incidence.  Using  Parseval's  theorem,  the  left 
side  of  (6)  can  be  evaluated  trivially  as 

J  E*(x,y) •  t!!,(x,y) dx  =  e'*"'” *•<” e'*'”*"  =  -k„ sin  0,..v)  ■  c  (8) 

One  usually  selects  basis  functions  that  have  closed  form  spectral  domain  representations;  therefore 
the  computational  effort  in  this  formulation  involves  the  numerical  evaluation  of  the  integrals  on 
the  right-hand  side  of  (6). 

3.  FREQUENCY-DOMAIN  IMPLEMENTATION  FOR  UWB  SIGNALS 
The  integrals  on  the  right  side  of  (6)  can  be  expressed  in  the  generic  form 

Imk  =  j  A(A;„ym,y*)e-<y*s-*Xe-'*-'^*dA:, 


(9) 
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where  A,  =  and  A,  =  x„-xl,.  Here,  y„  and  y*  locate  the  position  in  y  of  testing  function 

m  and  expansion  function  k,  respectively.  For  UWB  applications,  the  separation  /.  =  y  Aj+Aj 
between  expansion  and  testing  functions  will  range  from  zero  to  several  wavelengths.  Therefore, 
to  make  the  analysis  of  UWB  pulse  scattering  practicable,  special  considerations  are  required  for 
the  evaluation  of  integrals  of  the  form  in  (9). 


3.1.  Basis  functions 

Because  integrals  as  in  (9)  must  be  calculated  over  an  ultra-wide  bandwidth  and  the  efficiency 
of  such  integrations  determines  the  ultimate  speed  of  the  algorithm,  it  is  desirable  to  derive  a 
closed  form  asymptotic  expression  for  (9)  when  L  is  large  relative  to  wavelength.  As  will  be 
demonstrated  below,  it  is  possible  to  derive  such  an  expression  that  is  accurate  even  when  Z.  is  a 
small  fraction  of  a  wavelength. 

Success  in  this  endeavour  is  dictated  in  large  part  by  the  choice  of  basis  functions.  It  is  desirable 
to  use  basis  functions  with  simple  spectral-domain  representations.  For  this  reason,  the  complete 
domain  basis  functions  chosen  here,  which  do  not  explicitly  enforce  the  edge  condition,  are 


f,(x’)  =  sin 


'kv{x'  +  WI2)^ 

<w 

^ . 

•  i2 

(10) 


where  is  the  strip  width.  Note  that  in  (10),  the  explicit  y'  dependence  of  fi,(x'.y')  has  been 
suppressed.  This  y'  dependence  is  manifested  in  the  fact  that,  in  general,  the  strip  width  W  will 
be  different  for  each  strip  and  therefore  will  depend  on  the  layer  in  which  the  strip  is  located. 
The  spectral  domain  representation  of  /t(x')  is 

/*(k,)  =  jsin(k,W'/2)5*(/:.)  (11) 

for  k  even,  and 

A(k,)  =  cos(k,W/2)s*(kJ  (12) 

for  k  odd.  with  St(k,)  defined  as 

W 

The  spectral  representation  of  the  basis  function  therefore  consists  of  a  trigonometric  function 
which,  in  general,  varies  rapidly  with  respect  to  the  remaining  algebraic  expression  s,,.  It  should 
be  noted  that  the  commonly  used  triangular  subsectionai  basis  functions  also  can  be  written  as 
the  product  of  a  rapidly  varying  trigonometric  function  and  a  function  that  varies  slowly  in 
comparison.  Since  many  subsectional  basis  functions  are  usually  required  per  wavelength.^  they 
were  not  chosen  for  the  present  investigation  into  UWB  scattering. 


3.2.  Asymptotic  representation 

The  integrals  that  need  be  evaluated  can  be  grouped  into  two  types:  (i)  for  TM  waves 

=  rM  /A(*x)  )e“>v*^o-*i^e-i*.Adk,  (14) 

Jc 


and  (ii)  for  TE  waves 
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/n  =  WHO  f  yj=h(k,)  e  c  dk.  ( 15) 

JC  yfco  kl 

By  grouping  the  trigonometric  parts  of  the  basis  functions,  after  decomposition,  with  the 
exponentials  in  (14)  and  (15),  and  /™  each  can  be  expressed  as  a  sum  of  integrals  of  the 
form 


~  f  yjkl-kls„(,k,)s^(k,)e  dA,  (16) 

yvj«o)c 

and 

ATlf  =ja.Mof  (17) 

^kl-kl 

respectively,  with  n  either  0,  1,  or  -1. 

Tliese  integrals  can  be  approximated  to  a  high  degree  of  accuracy  using  standard  asymptotic 
techniques. ‘‘‘  It  is  convenient  to  introduce  the  following  change  of  variables; 


f,i  sin  0  =  .5,  + /ilV,  L,  cos  6  =  A...  L,  =  ^  =  k„sm  i  (18) 

which  transforms  (16)  and  (17)  into  the  following  simpler  equations  (see  Reference  14  for  the 
integration  path  in  the  4'plane): 


i 


=  T„.(Oi*(Ocos-4e-'«'«'«-‘'dC 


KT|  =  u»Po  UOi*(Oe 


(19) 

(20) 


with  n  =  koL,.  These  integrals  are  evaluated  most  efficiently  along  the  steepest  descent  path 
(SDP)  with  saddle  point  at  4  =  Q  '"*  By  performing  a  first-order  asymptotic  evaluation  of  (19) 
and  (2)  around  the  saddle  point,  one  finds 


K^-iDiiua  T.  ~s„(k,-- 


(21) 


KHI- 


2'Tr 

s„(k^=ko  sin  ©)  s^(  k,=k<,  sin  0 )  e'“  e'"  •* 


(22) 


Inspection  of  (21)  reveals  a  problem:  for  the  important  case  of  0  =  iTr/2  (expansion  and 
testing  function  on  the  same  plane),  (21)  predicts  =  0  for  ail  m  and  k.  This  is  because  an  x- 
directed  current  will  have  no  far-field  x-component  of  electric  field  along  its  axis.  Therefore  the 
first  order  approximation  in  (21)  is  only  valid  in  the  far  field.  A  second-order  asymptotic  evaluation 
is  in  general  quite  difficult.  However,  for  the  special  case  of  ©=±11/2,  a  second  order  approxi¬ 
mation  of  (19)  readily  yields 


Y ^^s„(A:,=A:oSin  0)s*(k^=knsin  0)  e  e'"  ''  (23) 

To  perform  higher-order  asymptotic  evaluations,  one  approximates  the  slowly  varying  part  of  the 
kernel  with  a  few  terms  of  its  Taylor  series,  and  therefore  one  must  differentiate  the  slowly  varying 
term  of  the  kernel.''*  This  is  trivial  for  the  kernel  in  (19)  when  4  =  if/2  because  the  derivatives 
of  5,„(4)  and  j*(4)  drop  out  in  view  of  the  vanishing  of  the  cos^(4)  term.  This  fact  was  used  to 
derive  the  approximation  in  (23)  for  4=»ir/2.  A  thorough  test  of  (21)-(23)  has  been  performed 
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for  a  wide  range  of  strip  widths  and  separations.  It  has  been  found  that  these  expressions  give 
good  agreement  with  the  numerical  evaluation  of  (14)  and  (15)  for  strip  separations  greater  than 
0-lXo-  As  an  example,  a  comparison  between  the  asymptotic  expressions  and  numerical  integration 
is  presented  in  Figure  2  for  Good  agreement  is  seen  over  the  entire  bandwidth.  Likewise, 
results  for  the  induced  currents  and  scattered  fields  computed  from  (21)-(23)  in  both  the  time 
and  frequency  domains  were  found  to  be  in  nearly  complete  agreement  with  results  calculated  by 
the  numerical  integration  of  ( 14)  and  ( 15)  (less  than  1  per  cent  difference).  As  a  practical  matter, 
it  should  be  pointed  out  that  one  must  use  L'Hopital's  rule  for  X,,  sin  S=‘±krrlW  to  get  good 
agreement  over  the  entire  bandwidth. 

In  summary,  use  of  the  asymptotic  expressions  to  replace  strictly  numerical  procedures  results 
in  a  tremendous  reduction  in  CPU  time  when  considering  UWB  plane  wave  scattering  from  a 
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Figure  2.  Real  and  imaginary  parts  of  I™  v.  where  k„  -  2Ttlk„  is  the  free  space  wavelength,  and  the  strip  separation 
equals  the  strip  width  W.  The  solid  line  represents  the  numerical  evaluation  of  (15)  while  the  dashed  line  represents  the 
asymptotic  approximation,  (a)  Real  part,  (b)  Imaginary  part 
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collection  of  strips.  This  is  true  especially  when  the  expansion  and  testing  functions  are  separated 
by  a  wavelength  or  more.  The  accuracy  of  the  asymptotics,  as  confirmed  by  comparison  between 
the  numerical  and  asymptoi.c  evaluations  of  the  integrals,  is  achieved  only  when  the  trigonometric 
parts  of  the  spectral  basis  functions  are  grouped  with  the  exponential  functions.  This  is  because 
the  trigonometric  parts  of  the  spectral  basis  functions  may  have  variations  equal  to.  or  greater 
than,  those  of  the  exponential  functions;  therefore,  they  must  be  combined  with  the  rapidly  varying 
parts  of  the  integrands  when  applying  asymptotics. 


3.3.  Integration  along  SDP 

For  expansion  and  testing  functions  separated  by  less  than  01\i,.  the  asymptotic  forms  in 
(21)-(23)  are  less  accurate  and  another  efficient  means  of  evaluating  ( 14)  and  ( 15)  is  required. 
For  such  cases,  ( 14)  and  ( 15)  are  evaluated  numerically  along  the  SDP. Although  this  is  obviously 
more  time-consuming  than  evaluation  of  the  closed  form  expressions  (21)-(23),  it  is  more  efficient 
than  performing  the  integral  along  the  real  k,  axis. 


3.4.  Self  term 

The  above-mentioned  techniques  are  useful  for  expansion  and  testing  functions  separated  in 
space,  and  hence  associated  with  different  strips.  For  the  self  terms,  however,  the  expansions  and 
testing  functions  occupy  the  same  strip  and  are  therefore  not  spatially  separated  (A,  =  dy  =  0). 
Again  ( 14)  and  (15)  can  be  reduced  to  a  sum  of  integrals  of  the  form  in  ( 16)  and  ( 17),  respectively. 
However,  since  \  =  0,  (21)-(23)  are  not  valid  for  n  =  (I  or  for  W  <  O  IX,,.  In  a  related 
problem,  it  has  been  demonstrated  that  integrals  of  the  form  (14)  and  (15)  vary  slowly  with 
frequency  when  A,  =  Ay  =  0.'^  This  can  be  unde,  stood  by  realizing  that  the  self  terms  sample 
essentially  the  near  fields  of  the  expansion  function,  and  the  near  fields  generally  vary  less  strongly 
with  frequency  than  their  far-field  counterparts. 

LWB  scattering  requires  very  fine  sampling  of  the  frequency  spectrum  in  order  to  furnish 
accurate  time-domain  results.  Realizing  that  (14)  and  (15)  vary  slowly  with  frequency  for 
A;r  =  Ay  =  0,  the  self  terms’  integrals  need  be  computed  only  at  points  along  a  relatively  coarse 
frequency  grid.  The  values  of  the  integral  between  points  along  the  coarse  grid  can  be  computed 
accurately  by  use  of  a  simple  extrapolation  procedure.'^  This  technique  of  computation  for  the  self 
terms  over  an  ultra-wide  bandwidth  has  been  applied  in  the  results  to  be  presented  subsequently, 
and  it  leads  to  significant  reduction  in  CPU  time. 


3.5.  Summary  of  integration  techniques 

In  summary,  for  expansion  and  testing  functions  separated  in  space  by  01\„  or  more,  the 
asymptotic  expressions  (21)-(23)  are  used  for  the  computation  of  the  reaction  integrals.  For 
expansion  and  testing  functions  separated  in  space  by  less  than  O-lX,).  numerical  integration  is 
performed  along  the  SDP.  The  expressions  (21)-(23)  require  virtually  no  CPU  time;  the  inte¬ 
gration  along  the  SDP  is  very  efficient  and  less  time-consuming  than  real  axis  integration.  The 
time-consuming  part  of  the  algorithm  for  UWB  applications  involves  the  computation  of  the  self 
terms.  Realizing,  however,  that  these  integrals  vary  slowly  with  frequency,  the  self  terms  are 
computed  only  at  points  along  a  relatively  coarse  frequency  grid.  All  self  term  integrals  at 
frequency  points  between  the  grid  points  are  efficiently  computed  using  linear  extrapolation.  It 
should  also  be  noted  that  due  to  reciprocity,  there  are  many  redundant  integrals  in  the  moment 
method  matrix.  Taking  advantage  of  this  redundancy  and  using  the  integration  techniques  summar¬ 
ized  above,  UWB  pulsed  plane  wave  scattering  from  a  collection  of  conducting  strips  becomes 
tractable.  It  is  believed  that  these  and  related  techniques  can  be  extended  to  other  classes  of  UWB 
scattering  problems. 
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4.  INVERSION  INTO  THE  TIME  DOMAIN 

In  order  to  obtain  accurate  results  for  UWB  pulsed  plane  wave  scattering,  it  is  necessary  that  a 
sufficiently  large  number  of  frequency-domain  points  be  used  before  Fourier  transforming  into 
the  time  domain.  The  integrals  which  vary  most  rapidly  with  frequency  will  be  those  associated 
with  the  most  widely  separated  expansion  and  testing  functions.  An  estimate  of  the  required 
number  of  frequency  points  required  can  therefore  be  found  by  examining  (21)-(23).  The  term 
in  these  expressions  having  the  most  rapid  variation  with  frequency  is  e''“.  Assume  that  E,  = 
is  the  largest  separation  encountered  in  the  problem  under  study,  and  that  is  the  highest 
frequency  component  needed  to  resolve  the  incident  pulse.  The  term  e  will  therefore  range 

from  1  to  where  is  the  free  space  wavenumber  at  If  N  frequency  points 

are  taken  per  period  of  oscillation,  then  Nk„  maxEn,a»/2rr  frequency  samples  are  required.  Here  is 
an  example  of  what  this  implies:  for  i-ma*  =  0-3  m.  =  28-3  rad/s  (100  GHz),  and  N  =  10. 
the  frequency  spectrum  must  be  discretized  from  0  to  100  GHz  in  100  MHz  increments.  This 
explains  why  it  is  essential  that  the  frequency-domain  results  be  computed  as  efficiently  as  possible. 

The  above  considerations  apply  only  to  time-domain  quantities  computed  directly  by  the  moment 
method  procedure  described  earlier:  the  time-dependent  currents.  To  compute  the  scattered  field, 
other  considerations  are  necessary.  In  the  far-held  approximation,  the  time-dependent  scattered 
fields  are  computed  from  integrals  of  the  form 

E(x.y.t)  =  jk(x.y.u))  c e’'"' db)  (24) 

where  r  is  the  distance  from  the  centre  of  a  given  strip  to  the  observation  point  (x,y).  The 
expression  h(x,y,o))  is  a  function  of  the  surface  currents  '.vhich  are  properly  described  in  the 
frequency  domain  by  the  discretization  procedure  discussed  above.  If  the  observation  distance  r 
from  a  given  strip  is  larger  than  (as  it  usually  will  be),  then  e"'*"'  will  vary  more  quickly 
than  e''*«^™«,  and  the  frequency  discretization  may  not  be  sufficient  to  describe  the  time- 
dependent  scattered  fields  accurately.  A  very  simple  procedure  can  be  used  to  overcome  this 
difficulty.  Equation  (24)  can  be  rewritten  as 


£,(j:.>’,7)  =  J  h{x,y.b))  e'“^  da>  (25) 

where  y  =  t-r/c.  The  expression  /i(x.y.a))  involves  only  the  basis  functions  and  the  Green's 
function  for  observation  of  the  fields  on  the  surface  of  the  strips.  The  discretization  required  to 
resolve  Eifx.y.y)  is  the  same  as  that  required  of  the  currents,  and  E(x.y,i)  can  be  found  easily 
by  shifting  Ei(x,y,y)  by  a  time  ric. 

It  should  be  noted  that  the  far-field  approximation  will  in  general  not  be  valid  for  significantly 
low-frequency  components  associated  with  a  given  incident  pulse.  However,  pulses  radiated  by 
practical  antennas  often  have  a  weakly  excited  low-frequency  spectrum  so  that  (25)  remains 
applicable  for  scattered  field  evaluation  in  most  observation  regions  of  interest,  including  those 
relatively  close  to  the  array  (but  'far'  from  each  strip).  All  time-dependent  fields  in  the  present 
study  have  been  computed  using  (25)  in  conjunction  with  the  FFT. 


5.  RESULTS 


5.1.  Frequency  domain 

To  the  authors’  best  knowledge,  there  are  no  results  in  the  literature  for  UWB  pulsed  scattering 
from  a  collection  of  strips.  Therefore,  to  check  the  accuracy  of  the  computer  code,  comparisons 
have  been  made  with  available  frequency-domain  results.  For  all  frequency  domain  calculations, 
a  total  of  12  basis  functions  was  u.sed  per  strip.  For  l.M  and  TE  incidence,  respectively,  Figures 
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3(a)  and  3(b)  show  the  normalized  currents  induced  on  three  eoplanar  strips  by  a  normally 
incident  plane  wave  (6,  =  0).  The  strip  widths  and  separations  are  X<,/4.  The  calculations  reveal 
good  agreement  with  data  computed  by  Cwik  and  Mittra/*  It  is  of  interest  to  examine  how  the 
trigonometric  basis  functions  resolve  the  edge  condition  for  the  case  of  TE  incidence.  From  Figure 
3(b)  it  is  seen  that  the  calculated  currents  oscillate  around  the  solution  found  when  the  edge 
condition  is  used  in  the  basis  functions”  (dashed  curve).  As  the  number  of  trigonometric  basis 
functions  is  increased  on  each  strip,  the  oscillations  become  more  closely  confined  around  the 
dashed  curve.  Although,  for  the  TE  case,  a  large  number  of  basis  functions  is  required  to  obtain 
adequate  convergence  for  the  currents,  it  has  been  determined  that  12  basis  functions  are  sufficient 
to  obtain  convergence  for  the  scattered  far  fields  (to  better  than  1  per  cent).  Therefore,  the  far 
fields  are  not  sensitive  to  the  above  discrepancies  in  the  surface  currents. 
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Figure  3.  Normalized  surface  currents  introduced  by  TM  and  TE  plane  waves  madent  vertically  on  three  sinps  of  strip 
width  and  separation  equal  to  Xo/4.  The  solid  line  represents  the  results  of  this  work  and  the  squares  represent  results 
from  Reference  4.  f  a )  TM  polarization  .(h)  TE  polarization .  The  dashed  line  was  computed  b  ;  including  the  edge  condition 

in  the  basis  functions*’ 
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Figures  4(a)  and  4{b)  show  the  scattered  far  field  due  to  a  TM  and  TE  plane  wave,  respectively, 
incident  at  0,  =  60“  upon  five  coplanar  strips.  The  strip  widths  are  O  IX,,  and  the  separation 
between  consecutive  strips  is  0  4X,).  These  results  can  be  compared  with  data  computed  recentiv 
by  Matsushima  and  Itakura  (Figures  6(a)  and  6(c)).'  It  is  difficult  to  transfer  the  data  accurately 
from  the  figures  in  Reference  5  (because  of  their  small  size),  but  excellent  agreement  is  noted 
upon  comparison.  It  is  worth  mentioning  that  the  basis  functions  in  Reference  5  for  TE  incidence 
satisfy  the  edge  condition  while  those  here  do  not.  Nevertheless,  the  agreement  between  Figure 
4(b)  and  the  results  in  Reference  5  is  excellent  for  all  observation  angles  except  very  near  i:'rr/2 
(abrupt  drop  in  the  pattern  of  Figure  4(b)).  The  discrepancy  over  this  very  small  range  of 
observation  angles,  which  may  in  fact  be  due  to  the  edge  condition  deficiency  of  the  basis  functions 
in  (10),  does  not  detract  from  the  utility  of  the  algorithm  employed  here. 


5.2.  Time  domain 

To  demonstrate  the  capabilities  of  the  numerical  code  for  UWB  pulsed  scattering  from  a  large 
collection  of  strips,  TE  and  TM  scattering  from  15  coplanar  strips  in  free  space  has  been  considered. 
The  strips  each  have  width  W  and  separation  2W/3.  The  incident  pulse,  normalized  to  the  width 
W,  and  its  frequency  spectrum  are  shown  in  Figure  5.  All  time-domain  results  are  plotted  as  a 
function  of  t.  where  t  is  the  time  required  by  a  plane  wave  to  travel  a  distance  W  in  free  space. 
Note  that  the  incident  pulse  has  a  temporal  length  shorter  than  t  so  that  it  is  capable  of  resolving 
individual  scattering  from  the  strip  edges. 

The  15  coplanar  strips  are  identified  as  follows;  the  centre  strip  is  defined  as  strip  0.  the  seven 
strips  to  the  right  of  the  centre  strip  are  labelled  1  to  ^  from  left  to  right,  and  the  seven  strips  to 
the  left  of  the  centre  strip  are  labelled  -1  to  -7  from  right  to  left.  The  r  =  D  time  reference  is 
defined  as  the  time  when  the  incident  pulse  reaches  the  centre  of  strip  0.  Figures  6(a)  and  6(bi 
show  the  scattered  fields  produced  by  TM  and  TE  plane  waves,  respectively,  incident  at  an  angle 
0,  =  20°.  The  scattered  fields  are  observed  at  a  distance  5()W/.3  directly  above  the  centre  of  strip 
0,  and  the  scattered  field  amplitude  is  normalized  to  the  peak  amplitude  of  the  incident  pulse.  In 
addition,  the  left  edge  of  each  strip  is  labelled  ’a’,  while  the  right  edge  is  labelled  'b'.  On  each 
strip,  the  incident  pulse  will  first  hit  edge  'b'  and  subsequently  edge  a’,  in  Figures  6(a)  and  6(b). 
the  travel  time  of  a  wavefront  along  a  straight  line  from  a  given  edge  to  the  observation  point  is 


1.0  .5  0  .5  10 


10  ,5  0  .5  10 

Figure  4.  Normalized  scattered  far  field  due  to  TM  and  TE  plane  waves  incident  on  five  coplanar  strips  of  width  (MKn 
and  separation  O  The  waves  are  incident  at  F),  =  Hf.  (a)  TM  polarization,  (h)  TE  polarization 
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Figure  5.  Incident  plane  wave  pulse  and  its  corresponding  normalized  frequency  spearum.  (a)  Pulse  shape  v.  tit.  where 
T  is  the  time  required  for  a  plane  wave  to  travel  a  distance  W  (strip  width)  in  free  space.  TTie  pulse  consists  of  three 
sections  of  a  sine  wave,  with  two  lobes  of  equal  amplitude  t  below  the  zero  axis  and  one  lobe  of  amplitude  2  above  the 
zero  axis,  (b)  Normalized  frequency  spectrum  v.  W/Xn 

labelled  for  all  strips  to  the  left  of  centre.  Note  that  as  one  moves  further  to  the  left  along  the 
strip  array,  the  scattered  pulses  from  individual  edges  become  more  distinct.  It  can  be  shown  that 
for  the  chosen  strip  distribution  and  incidence  angle,  the  scattered  wavefronts  from  the  edges  of 
strips  0-7  arrive  at  the  observer  at  nearly  the  same  time;  therefore,  these  signals  are  not  individually 
resolvable.  It  is  also  interesting  to  note  that  the  waveform  scattered  from  a  given  edge  is  different 
for  the  TE  and  TM  cases.  For  the  TM  case,  the  scattered  field  from  edge  ’b’  is  weaker  than  that 
from  edge  ‘a’,  whereas  the  opposite  is  true  for  the  TE  case. 

The  coptanar  strip  array  has  two  scales:  the  strip  width  and  the  strip  separation.  By  considering 
the  late-time  response,  and  focusing  on  the  time  delays  between  the  series  of  repetitive  waveforms, 
one  can  develop  a  scheme  by  which  these  two  scales  of  the  scar  .ring  cells  can  be  estimated  from 
the  data.  The  time  delay  between  consecutive  strong  and  weak  signals  gives  information  about 
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Figure  6.  Time-dependent  scattered  fields  due  to  TM  and  TE  plane  wave  pulse  in  Figure  S.  incident  at  6,  =  20°  upon  15 
coplanar  strips  of  equal  width  W  and  separation  2H'/3.  The  strips  are  labelled  as  follows:  the  centre  strip  is  strip  0.  the 
seven  strips  to  the  right  of  the  centre  strip  are  labelled  1  to  7  from  left  to  right,  and  the  seven  remaining  strips  are  labelled 
-I  to  -7  from  right  to  left.  The  left  edge  of  each  strip  is  labelled  'a',  the  nght  edge  is  labelled  ‘b'.  The  travel  times,  to 
the  observer,  of  wavefronts  from  the  edges  of  the  seven  strips  to  the  left  of  the  centre  strip  are  identified  by  arrows.  The 
time  reference  r  -  0  is  the  time  at  which  the  plane  wave  first  hits  the  centre  of  strip  0.  and  the  observation  point  is  SOW/3 
directly  above  the  centre  of  strip  0.  (a)  TM  polarization,  (b)  TE  polarization 


the  Strip  width,  while  the  time  delay  between  consecutive  large  (or  small)  pulses  gives  information 
about  the  strip  separation.  It  is  believed  that  the  insight  gained  from  such  simple  investigations 
will  be  useful  fc  understanding  the  scattering  of  UWB  pulses  from  a  more  general  class  of 
scattering  configurations.  A  detailed  analysis  and  explanation  of  these  and  other  results  obtained 
with  the  present  algorithm  will  be  submitted  separately  for  publication.'^ 

Finally,  some  observations  are  made  about  the  time-dependent  surface  currents  induced  on  the 
strips.  In  Figures  7{a)  and  7(b),  respectively,  are  shown  the  induced  currents  on  strip  0  for  the 
TM  and  TE  cases  investigated  in  Figure  6.  The  currents  are  plotted  as  a  function  of  time  at  three 
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Figim  7.  Time -dependent  currents  induced  on  the  centre  strip  of  the  IS-strip  anay  by  the  TM  and  TE  pulsed  plane  waves 
in  Figure  6.  Curves  I.  2.  and  3  are  for  positions  x  =  WI4.  e  -  0,  and  x  =  -WI4,  respectively,  on  a  strip  of  width  IF  with 
centre  at  x  =  0.  (a)  TM  polarization,  (b)  TE  polarization 


locations  along  strip  0:  at  x  ==  -WI4.  x  =  0.  x  =  W/4  (with  the  strip  centre  at  jr  =  0).  The 
incident  wave  induces  surface  currents  which  propagate  in  the  Jt-direction.  and  can  therefort 
expected  to  give  rise  to  resonances  between  the  strip  edges.  For  the  TE  case,  however,  the  indii 
currents  are  longitudinal  (z-direction)  and  are  therefore  expected  to  interact  less  strongly  between 
the  strip  edges.  These  expectations  are  confirmed  upon  examining  Figures  7(a)  and  7(b),  where 
the  late-time  oscillations  in  the  TM  case  are  much  more  pronounced  than  in  the  TE  case. 

The  15-strip  results  above  were  computed  on  an  IBM  6000  RISC  workstation  with  12  basis 
functions  used  per  strip.  The  time-dependent  data  required  calculations  at  2500  frequency  points 
(before  inversion  to  the  time  domain),  with  results  obtained  after  about  6-5  hours  of  CPU  time 
(for  an  average  of  less  than  10  CPU  seconds  per  frequency  point). 
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6.  CONCLUSIONS 

An  efficient  formulation  has  been  developed  for  the  analysis  of  UWB  pulsed  scattering  from  a 
large  collection  of  planar  strips  in  free  space.  By  the  spectral  domain  solution  strategy,  which  has 
been  summarized  in  section  3.5.  closed-form  asymptotic  approximations  have  been  derived  tor 
reaction  integrals  involving  expansion  and  testing  functions  separated  by  greater  than  01\«.  This, 
coupled  with  the  extrapolation  procedure  used  for  the  self  terms,  dramatically  reduces  CPU  time 
and  makes  .tie  analysis  of  UWB  scattering  tractable.  The  procedures  developed  in  this  paper  have 
been  applied  to  a  broad  parametric  study  of  UWB  scattering  by  different  arrangements  of 
strips,  to  various  processing  techniques  of  the  time-domain  data,  and  to  direct  and  quantitative 
interpretation  of  the  data  by  time-domain  wave  processes.  These  investigations  will  be  published 
separately.'^  It  is  also  intended  to  extend  the  algorithm  discussed  here  to  more  complicated 
environments  involving  dielectric  layers.  The  information  gathered  from  these  explorations  may 
find  use  in  the  interpretation  of  UWB  radar  data  from  large  periodic  and  quasi-periodic  multi¬ 
scale  environments,  such  as  ocean  waves. 
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SUMMARY 

We  show  the  derivation  of  parameters  in  the  frequency  domain  from  time-domain  data.  Far-held  character¬ 
istics  are  obtained  by  a  convolution  formula  with  the  harmonic  field  amplitudes,  which  are  obtained  via  a 
Fourier  transform  or  by  sampling.  The  electric  field  of  filter  ports  is  expanded  into  the  discrete  eigenmodes. 
By  this  method,  a  monochromatic  exact  open  boundary  can  be  formulated  and  the  fields  divided  into  the 
incident  and  the  reflected  part.  For  wide  band  operation  an  a  posteriori  error  correction  scheme  is  presented. 


INTRODUCTION 

The  analysis  of  electromagnetic  components  can  typically  be  subdivided  into  two  tasks.  First  the 
mathematical  problem  is  defined  and  solved  yielding  the  electromagnetic  fields  as  a  function  of 
one  temporal  and  three  spatial  co-ordinates.  The  second  task,  we  focus  on  in  this  paper,  consists 
of  reducing  and  filtering  the  result. 

One  common  method  of  eliminating  the  time-dependency  is  to  assume  harmonic  time-depen¬ 
dence.  In  the  case  of  constant,  time-invariant  materials  Maxwell’s  equations  are  decoupled  for 
different  frequencies  and  transform  to  quasistatic  differential  equations.  Furthermore  derived 
parameters  such  as  wave  amplitudes  are  used  to  describe  the  solution  in  order  to  obtain  a 
formulation,  which  is  analogous  to  a  discrete  network. 

The  direct  solution  of  the  frequency-domain  problem  using  finite  differences  or  similar  methods 
has  the  following  disadvantages.  The  equation  system  is  complex  and  therefore  twice  as  large  as 
in  the  time  domain.  When  calculating  near-resonances  the  algebraic  condition  may  become  very 
bad.  One  has  to  eliminate  the  spurious,  non-physical  solutions.  Also  we  have  to  repeat  the  solution 
process  for  each  frequency. 

The  alternative  is  to  calculate  the  time-domain  response  of  the  electromagnetic  component  and 
to  derive  the  frequency-domain  parameters.  In  the  following  we  show  the  calculation  of  far-held 
transforms  and  scattering  parameters  with  some  applications.  For  the  numerical  solution,  the  finite 
integration  algorithm  for  the  spatial  discretization  in  combination  with  a  leapfrog  scheme  for  the 
time  integration  was  used.*-^ 


FAR-HELD  CHARACTERISTICS 

One  disadvantage  of  finite  difference  and  finite  element  methods  is  that  all  computations  arc 
restricted  to  a  finite  grid.  Part  of  this  problem  can  be  overcome  by  introducing  radiation  boundary 
operators  simulating  an  infinite  mesh  size.  The  direct  calculation  of  far-held  characteristics  of,  for 
example,  antennas  is  still  infeasible,  since  the  grid  has  to  be  extended  to  distances,  where  the 
near  fields  have  ebbed  off. 

Therefore  we  have  to  strip  off  the  electromagnetic  fields  inside  the  grid  of  their  near-held  parts. 
The  far  field  then  can  be  written  as 


E,„(r.©,«b)  =  ~F(e,  <f>) 


(1) 


a  plane  wave  in  radial  direction  with  the  far-field  transform  F  (0.  <b)  as  the  directional  pattern. 
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Figure  1.  Integration  surface  around  radiator 

We  include  the  radiating  structure  in  a  volume  as  shown  in  Figure  I  and  substitute  the  radiating 
sources  by  equivalent  electric  and  magnetic  surface  currents,  given  by  the  tangential  electric  and 
magnetic  components  along  the  surface.  The  far  held  can  be  calculated  by  convoluting  the  surface 
currents  with  the  far  helds  of  elementary  eie''tric  and  magnetic  dipoles; ' 


F(0,<J)) 


4-11 


e,  X  (,  X  (n  X  B))  - 


1 

c 


X  (n  X  E) 


d.4 


(2) 


where  n  is  the  normal  to  the  integration  surface,  e,  the  normalized  radiation  vector  and  r'  the 
point  of  integration.  E  and  B  denote  complex  time  harmonic  electric  and  magnetic  field  amplitudes 
respectively. 

The  time  harmonic  field  amplitudes  can  be  obtained  either  by  using  a  time  harmonic  excitation 
and  sampling 


E  =  E(t„)-jE(f.,+ r/4)  (3) 

B  =  B(t„)  -  jB(f„  +  r/4)  (4) 

(T  denotes  the  length  of  the  harmonic  period  and  f„  has  to  be  a  time,  where  the  fields  have 
reached  their  harmonic  state),  or  by  an  on-line  Fourier  transform. 


FAR  FIELD  OF  A  CORRUGATED  HORN 

As  an  example  we  show  the  calculation  of  the  far  field  of  a  corrugated  horn.  The  structure  shown 
in  Figure  2  is  rotationally  symmetric  and  was  calculated  in  rz-geometry  using  65,000  mesh  points. 


Figure  2-  Geometry  of  the  corrugated  horn 
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As  excitation  a  pulsed  sine  wave  was  used.  The  geometrical  parameters,  normed  with  the  design 
frequency  ko  -  too/c.  are 


inner  diameter  of  waveguide 

bk„ 

3-25 

steepness 

d 

KT 

number  of  grooves 

N 

72 

depth  of  the  grooves 

r/X,, 

0-27 

distance  between  two  grooves 

01 

width  of  the  grooves 

UK 

0  067 

Figure  3  shows  the  time  response  of  the  radial  held  on  the  horn  axis.  The  time  harmonic  fields 
were  calculated  by  a  Fourier  transform. 

Figure  4  shows  the  far-field  transform  in  a  polar  plot  and  Figure  5  in  a  logarithmic  scale.  For 
comparison  the  results  calculated  by  R.  Erb^  by  a  mode-matching  technique  are  drawn  as  a  dashed 
line.  Both  results  show  good  agreement  except  near  90  degrees.  This  is  due  to  a  different  modelling. 
R.  Erb  assumed  radiation  into  a  halfspace  with  an  infinitely  conducting  screen,  whereas  here  a 
finite  structure  was  calculated. 


SCATTERING  PARAMETERS 

When  describing  multiports,  we  have  discrete  ports  and  a  discrete  spectrum  of  eigenmodes  which 
set  up  the  field  inside  the  waveguides  connecting  the  multiport  to  other  components.  So  we  use 
a  scattering  matrix 


[&i(ju>)b:(jw).,.6„(jw)|'  =  S(o,(jw)a2(ju))...a„(j(u))'  (5) 

to  describe  the  relationship  between  incoming  and  reflected  wave  amplitudes. 

The  transverse  electromagnetic  fields  are  described  in  the  frequency  domain  by  a  superposition 
of  incoming  and  reflected  waves 


Figure  3.  Field  in  the  horn  aperture  v.  time 
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Figure  5.  Logarithmic  plot  of  far  held  (dashed  line  results  by  R-  Erb) 


^  El' 

E(r,  <»)  =  Z  (aw( j<«>)e  +  fe„(jw)e^."“''l 


with  the  real  power  normalization  factor 


/■>(]“)  \'J  j  lE"  ^  H!{A:.y.ja))|dA 


(6) 


(7) 


For  the  further  derivation  we  use  the  weighted  wave  parameters  d^{ja»)  =  a^(jw)//„(ju>),  fe^(ju>)  = 
6w(jwV/»(j<*>>  and  formulate  the  time-domain  equation 
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E(r,  0  =  2  EXx,y,t)  •  [fl.(0  ’  +  bXi)  •  /»„(/, -2)1  (8) 


with  a  double  convolution.  The  frequency  dependencies  have  been  replaced  by  time  dependencies 
and 


P,(r.z)  =  ^J  e-^«H“’'ei’*-da)  (9) 

describes  the  waveguide  propagation.  In  homogeneously  filled  waveguides  we  can  find  frequency- 
independent  transverse  fields  E„(jr,y)  that  are  orthogonal. 

With  this  in  mind  we  look  at  the  following  model  of  a  transmission  line  (Figure  6).  z  =  0 
denotes  the  outermost  grid  line  at  the  waveguide  port,  z  =  8z  is  the  grid  line  one  step  inside  and 
A{t)  and  B(t)  are  the  transverse  electromagnetic  fields.  Both  can  be  expanded  in  terms  of  the 
discrete  two-dimensional  eigenmodes  EM))-  The  coefficients  are  composed  of  an  incident 

and  a  reflected  wave  amplitude: 

z  =  0 ;  =  d„(j<»)  -f-  S,.(  jui)  ( 10) 

z  =  8z:S„{jw)  =  fl„(ja>)P„(ja>)-'  +L(jw)£.(ja>)  (11) 


with 


P„(ju))  =  e-V)“'‘^  (12) 

We  can  realize  boundary  conditions,  that  are  exact  for  at  least  one  frequency  wm  when  we 
approximate  the  propagation  filter  Pw(j<^)  by  a  recursive  digital  filter  Pi(ju))  similar  to  classic  open 
boundaries'**  with  £„(j<i)M)  =  £»(jwM)-  With  this  filter  we  write  analogue  equations  for  the 
approximated  wave  amplitudes 

z  =  0:4,(j<ii)  =  d;(j<D)  +  fi'(juj)  (13) 

z  =  8z :  £Uj«)  = «:()«)  £:(jw)- '  +  §:(ju))  £:(ju,)  (h) 

and  formulate  a  recursion  for  the  unknown  reflected  amplitude 

5M)  =  PM)  *  (Sv(0  -  PM)  •  «:(01  (15) 

The  quality  of  the  boundary  condition  is  determined  by  the  filter  approximation,  but  we  can 
calculate  the  true  (weighted)  wave  amplitudes  d„(j(i>),  >n  a  correction  step  after  the  time- 

domain  field  calculation.  The  amplitudes  /4v(j(»),  Pv(j(»)  have  a  real  physical  meaning,  so  it  is 
possible  to  substitute  them  in  the  above  equations  to  yield  the  following  relationship 

z  =  0  z  =  Sz 


[ai(jw)g3(iu;)...o„(>w)j' 


A(t)  B(t) 

Figure  6.  Model  of  a  transmission  line 
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ra„(ju))i 

[  1  1  1 

~l 

[  1  1  ■ 

=  A(jw) 

.£v(jw)"‘ F:(ja)) 

.Svfjt*))] 

(16) 


By  this  technique  one  can  calculate  a  set  of  vectors  (ai(ja))  n^O***)-  S/.(j‘o))' 
b^(ja))..,b„(jui)]*  for  different  stimulations  (oKjuj)  oi(ju))...4' (jcu))'  to  solve  the  scattering  matnx 
5. 


BAND-STOP  FILTER 

As  an  example  the  results  calculated  for  a  waveguide  band-stop  filter  operating  in  the  X-band  are 
shown  in  Figure  7.  In  Figure  8  we  have  the  measured  and  in  Figure  9  the  calculated  transmission.' 
For  comparison,  results  calculated  by  a  mode-matching  technique  are  drawn  in  with  a  dashed 
line.** 

For  a  better  estimatior  of  the  error  the  curves  are  shown  in  Figure  10  using  a  range  of 
transmission  from  0  to  -  5  dB.  It  can  be  clearly  seen,  that  the  results  agree  within  0T5  dB. 


CONCLUSIONS 


In  this  paper  we  have  presented  two  methods  of  deriving  frequency-domain  results  from  time- 
domain  data.  Far-field  characteristics  are  calculated  by  applying  a  convolution  to  the  tangential 
time  harmonic  electromagnetic  components  on  the  surface  of  the  radiator.  These  fields  can  be 
either  sampled,  using  a  monochromatic  excitation,  or  obtained  via  a  Fourier  transform. 

For  the  calculation  of  scattering  parameters,  we  perform  a  mode  expansion  of  the  electromag¬ 
netic  fields  inside  the  ports.  The  i.iodal  coefficients  can  be  used  to  obtain  an  exact  open  boundary 
when  using  narrow-band  signals  and  contain  all  information  of  incident  and  reflected  wave 
amplitudes.  The  systematic  error  introduced  due  to  the  inexact  open  boundariej  in  a  wide  band 
range  can  be  compensated  by  an  a  posteriori  error-correction  scheme. 


Figure  9.  calculated  (ransmijaion  (dashed  line  results  obtained  by  mode  matching) 
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Figure  10  Calculated  transmission  (dashed  line  results  obtained  by  mode  matching! 
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SUMMARY 

The  Hilbert  space  representation  of  the  TLM  method  for  time-domain  computation  of  electromagnetic  fields 
and  the  algebraic  computation  of  the  discrete  Green  s  function  are  investigated.  The  complete  field  state  is 
represented  by  a  Hilbert  space  vector.  The  space  and  time  evolution  of  the  field  state  vector  is  governed  by 
operator  equations  in  Hilbert  space.  The  discrete  Green's  functions  may  be  represented  by  a  Neumann  series 
in  space-  and  time-shift  operators.  The  Hilbert  space  representation  allows  the  description  of  the  geometric 
structures  by  projection  operators,  too.  The  system  of  difference  equaiioi  governing  the  time  evolution  of 
the  electromagnetic  field  in  configuration  space  is  derived  from  the  operator  equation  for  the  field  state 
vector  in  the  Hilbert  space. 


1.  INTRODUCTION 

The  TLM  (transmission  line  matrix)  method  developed  and  first  published  in  1971  by  Johns  and 
Beurle  is  a  discrete  time-domain  method  for  electromagnetic  field  computation.'  '  In  this  paper, 
the  Hilbert  space  representation  of  the  TLM  method  is  presented  and  applied  to  the  algebraic 
computation  of  discrete  Green  s  functions.  The  Hilbert  space  representation  is  a  very  general  and 
powerful  concept  in  field  theory.^  Whereas  in  the  electromagnetic  theory  Hilbert  space  methods 
are  mainly  used  for  solving  the  field  equations,  as.  for  example,  in  the  moment  method.’  in 
quantum  theory,  the  fundamental  theoretical  concepts  have  been  formulated  in  Hilbert  space.*  ’ 

The  state  of  a  discretized  field  can  be  represented  by  a  vector  in  the  Hilbe.t  space.  The 
specification  of  the  mesh  node  connections  and  the  boundary  conditions  is  done  by  operators  in 
the  Hilbert  space.  The  Hilbert  space  representation  also  allows  the  description  of  geometric 
structures  by  projection  operators.  The  space  and  time  evolution  of  the  field  state  vector  is 
governed  by  operator  equations. 

In  field  theory,  field  propagation  in  spatial  domains  may  be  treated  using  Green's  functions." 
The  concept  of  Green's  functions  may  also  be  applied  to  discrete  time-domain  field  computation.'^ 
Discrete  time-domain  Green  s  functions  allow  the  modelling  of  the  relation  between  the  field 
values  on  the  boundaries  if  knowledge  of  the  field  in  the  spatial  domains  beyond  the  boundaries 
is  not  required. 

In  this  paper,  the  algebraic  computation  of  the  discrete  Green's  function  is  investigated.  Our 
approach  is  based  on  a  Hilbert  space  representation  of  the  space-  and  time-discretized  electromag¬ 
netic  field.  The  discrete  Green's  functions  may  be  represented  by  a  Neumann  series  in  space-  and 
time-shift  operators.  The  system  of  difference  equations  governing  the  time  evolution  of  the 
electromagnetic  field  in  configuration  space  is  derived  from  the  operator  equation  for  the  field 
state  vector  in  the  Hilbert  space.  First  results  are  presented  for  the  two-dimensional  case. 


2.  THE  TWO-DIMENSIONAL  TLM  METHOD 

The  electromagnetic  field  is  discretized  within  space  and  time.  The  space  is  modelled  by  a  mesh 
of  transmission  lines  connecting  the  sample  points  in  space.  The  field  computation  algorithm 
consists  of  two  steps: 
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•  The  propagation  of  wave  pulses  from  the  mesh  nodes  to  the  neighbouring  nodes. 

•  The  scattering  of  the  wave  pulses  in  the  mesh  nodes. 

In  the  following,  we  restrict  our  considerations  to  the  two-dimensional  case  with  the  transverse 
electric  field.  In  the  shunt  TLM  model,  voltage  wave  amplitudes  are  used  instead  of  total  voltage 
and  current.  The  voltage  wave  amplitudes  of  the  incident  and  the  reflected  waves  are  given  by 
and  The  left  index.  A.  denotes  the  discrete  time  co-ordinate  and  the  right  indices. 

m  and  «.  denote  the  two  discrete  space  co-ordinates.  We  consider  the  TLM  mesh  to  be  compiosed 
by  elementary  TLM  shunt  node  four-ports  as  shown  in  Figure  I .  where  each  of  the  four  arms  is 
of  length  A//2.  The  scattering  in  this  elementary  four-port  is  connected  with  the  time  delay  At. 
The  scattering  of  the  wave  pulses  is  described  by 


b, 

fi; 

=  S 

h. 

ttl.ti  A 

iia 

with  the  scattering  matrix  S  given  by 


(1) 


C) 


With  the  scattering,  a  time  delay  of  It  is  associated  and  therefore,  the  time  index,  k.  is  incremented 
by  one.  The  scattered  pulses  are  the  incident  pulses  of  the  neighbouring  elementary  cell.  This  is 
described  by 


l.^’l  <>i  •  I  /. 
rri  n-l 

„  .  I 


(3) 


3.  THE  DISCRETE  FIELD  STATE  SPACE 

In  the  TLM  model,  the  field  state  at  a  given  discrete  time  is  desenbed  completely  by  specifying 
the  amplitudes  of  the  four  wave  pulses  incident  to  each  mesh  node.  The  space  of  the  voltage  wave 
amplitudes  of  the  incident  and  the  reflected  waves  *0,  and  „  is  the  four-dimensional  real 
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vector  space  In  order  to  develop  our  formalism  in  a  more  general  way  we  introduce  the  four¬ 
dimensional  complex  vector  space  '€'*  for  representing  the  wave  amplitudes  *8^  „  and  *b„ 

In  order  to  describe  the  whole  mesh  state,  we  introduce  the  Hilbert  space  which  allows  the 
mapping  of  each  mesh  node  onto  an  orthonormal  set  of  base  vectors  of  The  time  states  are 
represented  by  the  Hilbert  space  5f,.  With  each  pair  of  discrete  spatial  co  ordinates  (m.n)  a  basis 
vector  of  is  associated  and  with  each  k,  a  basis  vector  of  is  associated.  We  now  introduce 
the  state  space  'X  given  by  the  Cartesian  product  of  ■€‘‘,  X„  and  X,. 

X  =  (4) 

The  space  X  is  a  Hilbert  space,  too.  The  complete  time  evolution  of  the  field  state  within  the 
whole  three-dimensional  space-time  may  now  be  represented  by  a  single  vector  in  X.  Using  the 
bra-ket  notation  introduced  by  Dirac."  the  orthonormal  basis  vectors  of  'X  are  given  by  the  bra- 
vectors  \k-,m,n).  The  ket-vector  {k:m.n\  is  the  Hermitian  conjugate  of  \k:ni.n).  The  orthogoi  ality 
relations  are  given  by 


{kx-.m,M,\k2-,m2.n2}  =  n’ 


(5) 


The  incident  and  reflected  voltage  waves  are  represented  by 


k>=  S  S  2 

rpfs-v  ft  -  • 


a\ 

Ui 

(lx 


‘k.m.n) 


k 


m.n 


and 


\b)=  f  2  2 

k,  -  V.  nt  -  -  n  ■  ■■  » 


b. 


bi 


\k-.m.n) 


(6) 


(7) 


in  the  Hilbert  space  M.  We  define  the  shift  operators  X.  Y  and  their  Hermitian  conjugates  X' 
and  Y’  by 


X\k:m.n)  =  iA.m-t- l./i> 

X*|A,/n.n>  =  \k:m- 1  ./i> 

Y|A;/n,n>  =  \k\m.n-^  1> 

\*\k-,m.n)  =  |A;m./i- 1>  (8) 

The  operators  X  and  Y  shift  the  field  state  by  one  interval  A/  in  the  positive  m-  and  /(-direction, 
respectively.  Their  Hermitian  conjugates  X*  and  Y’  shift  the  field  state  in  the  opposite  direction. 

We  define  the  time  shift  operator  T.  The  time  shift  operator  increments  A  by  1.  i.e.  it  shifts 
the  field  state  by  A/  in  the  positive  time  direction.  If  the  time  shift  operator  is  applied  to  a  vector 
\k-,m,n),  we  obtain 


T|A;w./i>  =  iA+  1  .m.n) 

We  introduce  the  connection  operator  P  given  by 

0  X  0  0 
X^  0  0  0 

r  = 

0  0  0  Y 
0  0  Y"  0 

With  the  connection  operator  P,  equation  (3)  yields  the  operator  equation 


(9) 


(10) 
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l^>>  =  r|a)  (11) 

describing  the  mesh  connections.  The  operator  F  is  Hermitian  and  unitary: 

r  =  r  =  r''  dZ) 

Therefore  we  obtain  from  equations  (11)  and  (12) 

k>  =  r  if>>  (13) 

We  now  express  equation  (1)  in  the  Hilbert  space  notation  by 

|f>>  =  TS|a>  (14) 

This  equation  describes  the  simultaneous  scattering  within  all  the  mesh  node  four-piorts  according 
to  Figure  1.  The  scattering  by  a  mesh  node  causes  the  unit  time  delay  St. 

Figure  2  shows  an  example  of  a  spatial  domain  within  a  TLM  mesh.  This  spatial  domain  is 
specified  by  a  given  set  of  mesh  four-ports.  A  spatial  domain  D  in  our  TLM  mesh  may  be  specified 
by  projection  operators.  We  define  the  domain  projection  operator  P„  which  projects  a  state 
vector  ia>  on  the  domain  D: 


P„ia>  =  !«>„ 

(15) 

This  projection  operator  may  be  written  in  dyadic  notation  as  the  sum  of  the  projection  operators 
on  the  nodes  belonging  to  the  domain  D: 

P/>  =  2  S  \l(-tn.n){k:m.n\ 

m&O  net) 

(16) 

In  the  same  way.  we  define  the  inner  domain  projection  operator  P, 
operator  by 

and  the  boundary  projection 

Pi  l«)  =  t«>i 

(17) 

P„  |a>  =  !<?>b 

(18) 

with 

P.  =  P,  P„ 

(19) 

Figure  2.  A  spatial  domain  within  the  TLM  mesh 


HILBERT  SPACE  FORMULATION  OF  TLM 


33 


Pb  =  PbPd  (20) 

Pb  +  Pi  =  Po  (21) 

The  inner  domain  projection  operator  projects  the  circuit  space  IK  on  the  inner  ports  of  the 
domain  D  (Figure  3).  Since  the  projection  operator  P,  and  the  connection  operator  F  arc 
commuting,  i.e. 


|Pi  r]  =  0  (22) 

we  obtain 

l<>>i  =  r|a),  (23) 

Applying  diakoptics  to  TLM  structures  requires  the  computation  of  the  wave  pulses  scattered  at 
the  domain  boundaries.  The  initial  conditions  or  boundary  conditions  are  given  by  the  wave  pulses 
incident  on  the  boundary  ports.  We  apply  the  projection  operators  P,  and  Pr  in  order  to 
separate  the  field  states  {u)  and  lb)  into  the  inner  field  states  \a)i  and  16),  and  the  boundary  states 
\a)n  and  (6>b.  From  equation  (14)  we  obtain 


|(t)B  ~  T  Sbb  |u)b  +  T  Sm  ju), 

\h),  =  T  S,B  \a)n  +  T  S„  |u),  (24) 


with 


Sbb  ~  Pb  S  Pb 
Sbc  ~  Pb  S  P, 

Sm  =  P|SP„  (25) 

S„  =  P,  S  P, 

Using  equations  (23)  and  (24).  we  eliminate  the  inner  domain  states  |a),  and  16),  and  obtain 

|i>)B  =  lTSBB+TSB,(l  -  FTS,,)  'rTS,„llfl)B  (26) 

This  is  the  relation  between  the  incident  and  .scattered  boundary  state.  It  describes  the  evolution 
of  the  boundary  field  state  without  knowledge  of  the  inner-field  state.  It  has  to  be  considered  that 
the  operator  equation  (26)  is  non-local  with  respect  to  both  space  and  time.  We  expand  the 
operator  (l-F  T  S,,)"'  into  a  Neumann  series'"  "  and  obtain 


Figure  3.  The  inner  pons  of  a  TLM  domain 
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(i-rTs„)-'  =  VT'(r  s,,)' 

/■•(j 


(27) 


Inserting  this  into  equation  (26)  yields  the  boundary  state  evolution  equation 

\h>B  =  G|fl>B  (28) 

with  the  boundary  held  evolution  operator  G  given  by 


G  = 


tSbb  +  Sb,(  iT'’^(rs,.)')rs„ 


(29) 


The  boundary  held  operator  G  gives  the  relation  between  the  boundary  state  vector  jals  rep¬ 
resenting  the  wave  pulses  incident  on  the  boundary  and  the  boundary  state  vector  {6)0  representing 
the  wave  pulses  reflected  through  the  boundary.  Equation  (28)  is  the  general  formulation  of  the 
boundary  element  problem  in  the  Hilbert  space.  Since  the  Neumann  series  is  an  inhnite  geometrical 
series  in  space-  and  time-shift  operators,  the  boundary  held  operator  is  non-local  with  respect  to 
space  and  time. 


4.  THE  DISCRETE  TWO-DIMENSIONAL  GREEN  S  FUNCTION 

As  an  example,  we  derive  the  discrete  Green  s  function  for  the  half  plane.  The  discrete  Green  s 
function  for  the  half-plane  is  given  by  the  projection  of  the  boundary  state  evolution  operator 
equation  (28)  onto  conhguration  space  for  a  point-like  initial  state  la)B.  The  half-plane  (Figure  4) 
is  dehned  by  the  domain  projection  operator  Pf,  given  by 


Pw  =  S  2  (30) 

k.n  m  all 

As  in  the  shunt  TLM-model,  voltage  wave  amplitudes  instead  of  total  voltages  are  used,  a  new 
Green's  function  for  wave  amplitudes  has  to  be  dehned.  For  a  boundary  problem,  the  discrete 
Green's  function  is  dehned  by  the  convolution 


Figure  4.  The  homogeneous  iwo-dimensionai  half-space 
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k^n  A  -  A'^n  ~  n'  k'^n' 

k'  .ft' 

=  kG„  •  (31) 

with 

neB  (32) 

where 

.../iv}  (33) 

denotes  the  set  of  N  boundary  nodes. 

*fl«  is  the  column  vector  of  the  incident  impulse  functions  at  the  time  k^t.  i.b„  is  the  column 
vector  of  the  scattered  output  wave  pulses  at  the  time  kSt.  tG„  is  the  discrete  Green  s  function 
for  an  arbitrary  boundary  with  A/  boundary  nodes.  It  describes  the  relation  between  the  incident 
and  the  scattered  wave  amplitudes  in  the  boundary  ports. 

For  the  half-plane,  the  boundary  is  given  by  m=0  and  «  =  l.().l .. ..‘.c.  Therefore  equation 

(31)  yields 

A,=  ^  *  ,  (34) 

The  boundary  state  evolution  equation  (28)  may  be  expressed  by  the  discrete  Green  s  function, 
equation  (34),  via 

|/>>o  =  G|a)„  (35) 

where  the  boundary  field  evolution  operator  is  given  by 
10  0  0 

G-  y  ”  i;  i  ,.uG„-„  \k-Q.n){k-:0.n'\  (36) 

0  0  0  0 

In  order  to  calculate  the  Green  s  function  for  the  boundary  of  the  half-plane,  we  start  from  an 
impulsive  excitation  at  /i'=0,  ^'=0  given  by 

I 

I  10;0.0)  (37) 

0 

and  obtain 

1 

|ft)n=  i  i  *GJ^;0,n)  (38) 

,i  —  *.  4-  ^ 

0 

Our  result  will  be  the  Green  s  function  *G„. 

Mapping  equations  (28)  with  (29)  and  (37),  (35)  with  (36)  and  (37)  to  configuration  space  by 
multiplying  both  equations  from  the  left  side  with  (A:;0,«|,  we  obtain 
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{A:;0,wj6)B  —  ^0,n|SBt  (T  Sn)*  ■  F  S|h 


(k-,0.n\b)n  = 


where  we  have  used 


(AIT'IO)  =  8,., 


so  that  only  one  term  of  the  Neumann  series  in  equation  (29)  contributes  to  jC,,.  if  we  restrict 
ourselves  to  ks=2. 

With  this  result,  we  consider  equation  (29)  and  formulate  the  main  part  of  the  problem  that 
means  the  operator  (T  f  Sn)'  recursively 

!«>».,  =Trs„|«),  (42) 


[o)k  -  2#  ” 


With  the  projection  operators  Pi  and  Pu  given  by 


0  ()  0  (•’ 
0  10  0! 
0  0  I  Oi 
0  0  0  1  : 
10  0  0 
0  10  0 
0  0  I  0 
0  0  0  1 


2  \ki).n){kS).n\ 


2  S  \k.m.nHk:mju 


1  0  0  0 

0  0  0  0  „  , 

Y  \k-.i).n)<.k:0.n\ 

0  0  0  0  rr,' 

0  0  0  0 


it  yields  for  the  operators  Sbb.  Sbi,  S|b  and  Sn: 


-i  0  0  0 

0  0  0  0,,^, 

A  A;0,n)  (/(r;0,« 

0  0  0  0  ri 

,0  0  0  0 
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We  obtain 


and 


Defining 


Sb(  — 


Sib  ~ 


S„  = 


0  §  1  J 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
i  0  0  0 
i  0  0  0 
i  0  0  0 
'  -i  i  J 
i  -i  i 
1  i  -i 
.  i  ^  i 
0  0  0 
0  -i  -i 
0  i  -i 


2i*;0.n>(*;0.n| 


XM.rt></t;0./.| 


i 

i 

-i 


0  J 


0 
.! 

I 

*  -i 


2  S  |A’;w,n)  (A;m.n| 


2  |*;0./i)(/(:;0,ni 


'  1 

■  |1;1.0> 

|a)*..,=TrS,B 

0 

|0;0.0>  =  ] 

0 

0 

il;0.1) 

.  0 

|1:0.-I) 

1 

■  0  i  i  i 

0 

=  <*;0.nlTSB,  K  -,  = 

0  0  0  n 

0 

0  0  0  0 

.0, 

,0  0  0  0 

<<c;w,n|<i>*  = 


A 

L  J  m.n 

1 

0 

0 

1 

0 

0 

0 

,  0 

0 

0 

0 

+ 

1 

0 

0 

1 

(47) 


(48) 


(49) 


(50) 


(51) 


(52) 
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and  assigning  k=k'  —  \.  we  obtain  the  following  system  of  partial  difference  equations  by  mapping 


equation  (42)  to  configuration  space 

“  +  *(fl4]m./i)  (53) 

+  i  l./t  A[^4jm +■  I ./?)  (54) 

4  1  ^  i  (A{^l|fn.n  kt^2\m.n  *[^3]m.rt  Jk{^4)m.n)  (55) 

Ac  +  l(Ui],fi.n- 1  “  i  (A(^ljm.n  4-I^2jm,#l  “  *(03]^,,,  4[^4lm./t)  (56) 


for 

it  =  0.1.2... .* 
n  =  -»....-  1,0,1. ...X 
m  =  0,1,2.. ..X 

The  initial  conditions  are  given  by 

nf^iji.o  “  i  ~  5  (57) 

o[a4)().-i  =  1  all  other  w,n;  =  0  (58) 

As  the  space  is  not  bounded  with  respect  to  n.  we  only  need  boundary  conditions  concerning  m. 
One  boundary  value  is 


=  0  (59) 

As  we  have  a  system  of  second  order  concerning  m.  we  need  another  boundary  condition. 
Therefore  we  apply  the  Sommerteld  radiation  condition. 

From  equation  (51).  we  obtain  for  the  Green's  function  for  it  s'  <1 

“  *|Uo)<l.n  ~  1  (*(<*2]t).o  *[U3j«.„  +  *(i*4]l).n)  (60) 

This  system  of  partial  difference  equations  can  be  solved  by  transforming  it  to  frequency-  and 
momentum-space.  Concerning  n,  we  consider  ^[a,],„  „  as  the  Fourier  coefficients  of  the  function 

S{*["iU.«}  =  a'4-(C)]cn=  S  *kL..r"  (61) 


with 


^  =  exp2TrjA  (62) 

Concerning  m  and  k.  we  apply  the  Z-transformation'-* 


Z{  *[A, (?)!,„)  =(S,(4,v)j„=  2  *(A,(4)i„v  * 

<lc=0 


(63) 


with 


V  =  exp  2iTj/ 


(64) 


and  in  analogy  to  equation  (63)  the  Z-transformation  with  respect  to  m 
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Y{[S,(€.V)U}  =  (C{fn,v)|  =  t  - 


with 


=  exp  jA: 

As  we  have 

for  a  shift  in  the  positive/negativc  direction  of  n  and 

Z  {*..{A,(OU}  =  -  v„[4, (?))„, 

for  a  left-shift  of  k  and  m.  we  obtain 

2i>T)C|-V  =€,-€;  +  C,  -t-  Cj 
2;C,  -F  (Bol„  =  -  C,  +  C,  -F  C,  +  G 
2v^C,  -  V  =  C,  +  G  +  C,  -  G 
2'G-»'  =  C,  +  G-G  +  G 

where  we  have  used  the  abbreviation 

(®oj«  (^;ln  + 

This  system  of  algebraic  equations  has  the  so'ution 

C,  =  —  (2fv'  -  -  -qv-  -  2^v^  -  q^-v-  +  3£qv) 

+  ^-^5^  -  qv  -  q5=v  +  |q) 

^2  =  +  q*4^P^  -  -2|qV  +  ^qv) 


[flo)o 

JV 


(24qV  -  Iqv^  -  q-v^  -  q^^’v^  +  |q) 


Cj  -  ^  (2qV'  -  25qv^  -  qV  -  ^q-v^  -  qv=  -F  ^qi/) 


[^o)o 

/V 


(q  V  -  ^q^v  -  q  V  -F  4q) 


G  =  ^  (2?^qv^  -  V  -  2^ir>  +  ^q V  -  f  qv*  +  ^qv) 

_  (|2-q2v^  -  nf^v  -  ^q^v  -F  ^q) 

with  the  nominator  N  given  by 

N  =  2q^  (^v  —  -F  2q(24v''  —  +  v  -F  —  2^)  -F  2^w  —  2^v^ 


(65) 


(66) 


(67) 


(68) 

(69) 


(70) 

(71) 

(72) 

(73) 


(74) 


(75) 


(76) 


(77) 


(78) 


(79) 
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For  the  inverse  transformation  with  respect  to  m.  we  consider  the  function  C  4  and  rewrite  equation 
(78): 


C4 


0-  -  -T)  (  V  +  '  -  5  -  )  +  1 


I-  )2ev'-2v--£v-l 

*  *  L»  > 


2(l-v-) 


n-  -  2ti  V  +  ^  -  ! 


We  assign 


and 


with 


It  yields 


cosh  a  =  i’ 


1  £ 


1 

2£ 


sinh  a  =  e  \  cosh-  a  -  1 


1  for  JJ{a)  s:  {) 

-  1  for  iyi{o}  S'  0 


tl^-nc(«ha  |v  +  [fl,,],,)  (1  -  £v) 
-  —  271  cosh  a  +  1  2(1-1'-  ) 

Tt  ^inh  n 


with  the  function  T(i,v]  given  by 


(80) 


(81) 


(82) 


(83) 


(84) 


7'(£.v) 


-  iv-  +  i£-v'  -  i£v  ~  2 


2{l-v^)^(v+’-f-~J  -1 


+  l8o)<> 


£v-  a-  iv  +  i^=v  -  H  ~ 


1 


2(1  -V=)  s/(v+  '  -  iV-  1 


(85) 


The  following  correspondences  are  valid:' 


Y  {cosha/n} 


n*  ~  T)  cosh  a 


n’  —  27)  cosh  o  ■+• 


(86) 
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Ti  sinh  a 

Y  (sinhom)  =  -  r - . —  ( 

'  IT  -  2ti  cosh  a  f  I 

with  m  = 

Writing  out  the  sinh(a/n)-  and  cosh(am)-function  in  exponential  terms,  we  obtain  for  (Sjj 


ID)  in  I  . /(P  .u 

[flaU  =  b2  exp(  +am)  |  20 

^  ,n _ / 


!exp(-am) 


2(1  -  wO 


cm.v) 


The  Sommerfeid  radiation  condition  means  that  for  Fla}  5?  0,  the  terms  with  exp(-am)  and 
for  R|a}  ^  0,  the  terms  with  exp(+aw)  must  vanish  because  in  passive  media,  exponentially 
growing  solutions  do  not  correspond  to  physical  solutions.  In  both  cases  we  obtain 


0-  -r{ «,.]..)( 1  -  £>•) 
2(1  -  !•-) 


and  for  (fl,i|ii 


[B„{6,v)i„  = 


2v-(  v'  -  1 )  y  I  »■  •  -  cos  t)  I  - 

V-  -  2vcos  H  *  1 

2i-’  -  2v'  cos  (J  -  v'  ^  V 
V-  -  2v  cos  a  -  1 


for  m=l).1.2....*  with 


cos  6  =  cos  2  rr.V  =  i  £  + 


We  only  need  the  boundary  value  {S„(0,v)jii  to  calculate  i.G„.  because  with  equation  (60)  it  yields 

,,2C„  =  4a.,l„,„  =  IX-'  (Z-’  ;lS„(4.v))„})  (92) 

for  k  2^  0.  Of  course,  the  result  for  (Sn)m  the  same,  if  we  consider  the  functions  C,.  CT  or  C,. 
The  transformation  back  to  time-space’’  can  be  achieved  by 


where  C  is  a  closed  curve  in  the  complex  v-plane  which  surrounds  the  unity-circle.  Integration 
must  be  taken  in  anticlockwise  sense. 

Owing  to  the  orthogonality  relation 

•i 

(^(N))’’"’ tl/V  =  (94) 

.0 

the  inverse  transformation  from  momentum  space  to  configuration  space  concerning  n  is  given  by 
X-‘  {*1/1|(0)U}  =  LkU  =  2!^  J^"  JT,(e)Uexp(j0n)de  (95) 


As  BofC.*')  's  function  of  0,  it  yields  with  equations  (93)  and  (95)  for  i,(a,i)n., 
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*(^lo.n  =  I  -  [  (s<.(e.p)]..  COS  «e  p*  “ '  de  dv 

in/  )c  IT 


3i/i (u^  -  \)  +  y  -  cos e j  -  1 


ir  —  2v  cos  P  +  1 


=  J-Iif 

2''fj  Jc-^Jo 

“  I  ~  T1  ne  V*  -  ‘  de  dp 

2irj  Jf  ir  J„  P"  -  2p  cos  9  +  1 


■  cos  r0  V*  ‘  d6  dp 


(96) 


where  we  have  used 


j  (analytical  function)  dp  =  0 

For  the  evaluation  of  the  two  integrals  /|  and  /,.  we  can  restrict  ourselves  to  the  case  n  »  0 
because 


With'-* 


for  p-  >  1.  «  35  0  and 


we  calculate 


cos  nO  de 
p-  -  2pcos  e  +  1 


up 

p- 


I 

2-irj 


dp  =  6,,„ 


^2  ~  ~  8*  n-l 


For  the  integral  /,,  we  apply'^ 


=  i  p'p,(x) 

y  1  +  2pjr  +  p- 


for  IpI  <  1  and  |?|  «  1  and'* 


i-p- 

1  -  2pcosjr  +  p~ 


X 


=  1+22  ^ 

i-l 


(97) 


(98) 


(99) 


(100) 


(101) 


(102) 


for  IpI  <  1. 

With  the  integral  representation  of  the  Legendre  polynomials'^ 


1  ^  r  df 

2-^)  jc-  ^{2  -  2xt  +  1 

where  C  is  a  closed  curve  in  the  complex  p-plane  which  encircles  the  unity-circle  in  anticlockwise 
sense,  we  obtain 


Pn(x) 

0 


forn  s  0 
■or  n  <  0 


(103) 
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l\  ~  k-\ln 

k+l 

;-0 

k 


(104) 


where  the  function  */„  is  defined  by 


P,(cos^B/2)  /’t-,(-sin‘ e/2)cosn0d0 


(105) 


To  calculate  */„,  we  expand  iP,(cos"  6/2)  and  P*-/(-sin^  0/2)  in  terms  of  cos-  6/2  with  the  help 
of  References  18  and  19: 


/’*_,(-sin-0/2)  =  (-1)*  '  X  ^  ^)  (-))'■  (cos- 6/2)-' 

1  /7\  /^f  - 

P,(cos‘  0/2)  =  2,  S  (-  I)'  (^)  (“  /  )  (cos=  6/2)='-' 

We  substitute  6  =  2>',  apply  the  integral^'’ 

2="“  U-m/ 


(106) 

(107) 


‘nl 

Jo 


cos^'’ t  cos  2mt  dt  = 


for  n  3  m 
forn  <  m 


(108) 


and  obtain  for  the  function  */„; 

*  (//:|  *  - 1 


.'--22  2 

/=0  ,«0  r-o  '  *  ’ 

Combining  the  two  integrals  yields  an  algebraic  expression  for  the  discrete  Green’s  function 


k(jn  ~  ^^k.n+l  i  i  H +  l7n  i*-j4i 


i-0 

k-2 

2  ^  A— 2-/'7/f  +  I  — /  —  1  k-2-'fJn-  l—f  ^  k~2~fJn~3~f 
i-0 


(110) 


for  /i=0,l,2,...’c  and  ^ -2.3,4.. ..^c.  Because  of  equation  (97)  we  have  for  n  «  0; 

kG.„  =  ,G„  (111) 


As  already  remarked,  the  general  Green’s  function  for  an  excitation  at  the  time  k'  in  the  boundary 
node  n'  is  obtained  by  the  transition 


kG„-*  k-k  G„-„- 

In  Figure  5,  *G„  is  depicted  for  n=-9,  ...,-1,0,1 . 9;  ^=1.2,  ....  10. 


(112) 
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.2 


.0 


Figure  5.  Vjiucs  ot  ihe  Green  s  function 


5  CONCLUSIONS 

In  this  paper,  the  TLM-method  has  been  represented  in  the  Hilbert  space.  The  Hilbert  space 
formulation  allows  the  derivation  of  general  algebraic  expressions  for  the  held  evolution  without 
having  regard  to  the  individual  geometric  conditions.  Geometrical  structures  may  be  described  in 
a  general  way  by  projection  operators.  A  further  advantage  of  the  Hilbert  space  formulation  is 
that  the  powerful  methods  of  functional  analysis-'  can  be  applied.  General  investigations  of 
operator  equations  concerning  the  existence  and  the  convergence  of  the  solutions  are  simplified. 

The  Hilbert  space  formulation  was  used  for  deriving  general  algebraic  expressions  for  the  space 
and  time  discrete  field  evolution.  From  the  local  operator  equations  governing  the  time  evolution 
of  the  field  state  vector,  the  non-local  operator  equation  describing  the  time  evolution  of  the 
boundary  state  vector  was  derived. 

First  applications  of  this  method  were  demonstrated  in  calculating  the  discrete  Green  s  function 
for  the  half-plane. 
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SUMMARY 

In  this  paper  a  new  TLM  model  for  the  two-dimensional  wave  equation  is  introduced.  It  is  synthesized 
directly  from  a  FD  algorithm.  The  FD  algorithm  is  second-order-accurate  in  both  space  and  time,  and  is 
explicitly  time-stepped.  The  spatial  derivatives  in  the  FD  algorithm  are  approximated  by  the  weighted 
combination  of  two  standard  central  difference  stencils,  one  oriented  as  usual,  the  other  rotated  by  45°  with 
its  arms  extended  by  a  factor  of  (2)*'^.  The  TLM  model  is  realized  as  the  weighted  connection  of  two  original 
models  (with  the  same  geometrical  configuration  as  the  FD  algorithm).  The  weighting  in  the  TLM  model  is 
accomplished  by  using  a  variable  intrinsic  impedance  for  specific  elemental  transmission  lines.  The  FD  and 
TLM  methods  possess  identical  dispersion  relations  if  the  former  is  operated  at  its  upper  limit  of  stability. 
Therefore,  under  these  conditions  both  represent  identical  models  for  the  simulation  of  wave  propagation. 
The  propagation  characteristics  of  the  new  model  are  investigated  and  the  conditions  for  approximate 
numerical  isotropy  are  provided.  The  numerical  implementation  (scattering  matrix  and  transfer  event)  is 
described.  To  validate  the  new  model,  the  calculation  of  cutoff  frequencies  of  various  modes  in  rectangular 
waveguide  is  performed.  Comparison  with  analytical  results  (for  an  unfilled  waveguide)  and  other  numerical 
results  (for  a  waveguide  partially  filled  with  a  dielectric)  validate  the  implementation  of  the  model. 


1.  INTRODUCTION 


The  numerical  techniques  discussed  in  this  paper  are  capable  of  solving  arbitrary  two-dimensional 
electromagnetic  field  problems.  If  problems  independent  of  the  ;-direction  are  considered. 
Maxwell’s  equations  are  reduced  to  two  independent  sets,  one  of  which  is  given  by. 


da:  ^  r)l 


dy 

d//,  d//, 

dx  dy 


(t£-  +  c 


aE, 

dr 


(la) 

(lb) 

(lc) 


where  Ep  and  Hp  are  the  electric  and  magnetic  fields,  respectively  (with  p  =  x,  y,  or  z)  and  c, 
p.,  and  a  are  the  permittivity,  permeability  and  conductivity  of  the  medium  of  interest,  respectively. 
Equations  1  can  be  combined  to  yield  the  two-dimensional  wave  equation  in 


dx^  ^  dy^ 


op 


d£, 

df" 


+  ep 


d-£, 

df= 


(2) 


The  numerical  techniques  presented  in  this  paper  are  developed  from  discrete  approximations  to 
(2)  rather  than  (1). 

Johns  and  Beurle  introduced  the  transmission-line  matrix  (TLM)  method  in  1971  as  a  technique 
which  utilizes  the  equivalence  of  voltages  and  currents  on  transmission  lines  to  electric  and 
magnetic  fields  in  space.'  An  orthogonal  grid  of  transmission  lines  represents  a  physical  model 
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which  approximates  (1)  or  (2).  Hoeter*^  presented  recent  applications  and  extensions  of  the 
method.  Another  numerical  technique  used  extensively  in  the  computational  electromagnetics 
community  is  the  finite-difference  time-domain  (FD-TD)  method  introduced  by  Yee^  and 
extended  by  others.’  Both  the  TLM  and  FD-TD  methods  are  capable  of  providing  approximate 
solutions  to  the  time-dependent  form  of  Maxwell  s  equations.  In  their  most  basic  forms,  both 
utilize  regular  rectangular  grid  structures  and  explicit  time-stepping.  Under  certain  circumstances 
both  methods  represent  identical  models  for  wave  propagation.  For  all  cases  in  which  an  equival¬ 
ence  between  a  TLM  model  and  FD  algorithm  has  been  established,  the  TLM  model  corresjxmds 
to  the  FD  algorithm  when  the  latter  is  operated  at  a  specific  location  in  its  stability  range.®  '* 

Recently,'*  the  equivalence  of  the  original  TLM  model'  and  the  two-dimensional  Yee  algorithm'* 
is  established.  In  Reference  7  the  equivalence  of  the  three-dimensional  expanded  node'*  and  the 
three-dimensional  Yee  algorithm*  is  established.  In  Reference  8.  models  of  (2)  based  on  hexagonal 
(rather  than  rectangular)  computational  grids  are  investigated,  a  TLM  model  is  presented  and  its 
equivalent  FD  algorithm  derived. 

In  general,  the  finite  difference  (FD)  method  can  be  applied  in  various  va'S  to  approximate 
(2).  Grid  structures  and  the  accuracy  of  the  difference  tormulas  can  be  varied,  and  different  time¬ 
stepping  schemes  can  be  used.  The  purpose  of  this  paper  is  to  synthesize  an  equivalent  TLM 
model  directly  from  a  FD  approximation  of  (2).  The  general  approach  can  be  extended  to  the 
synthesis  of  other  TLM  models  from  FD  algorithms. 

In  the  following  section,  the  FD  algorithm  is  presented  us  a  weighted  connection  of  two  Yee 
algorithms, •*  one  oriented  as  usual  (arms  of  the  spatial  stencil  oriented  along  the  x-y  axis),  the 
other  rotated  by  45®  with  its  arms  extended  by  a  factor  of  (2)'  ^  The  dispersive  characteristics 
and  stability  criterion  of  the  algorithm  are  derived.  In  section  the  equivalent  TLM  model  is 
presented.  Based  on  the  relationship  established  in  Reference  6.  the  equivalent  TLM  model  is 
constructed  from  an  interconnection  of  two  original  models.  One  oriented  as  usual  (elemental 
transmission  lines  oriented  along  the  x-y  axis),  the  other  rotated  by  45®  with  its  arms  extended 
by  a  factor  of  (2)'^^.  The  weighting  is  accomplished  through  the  use  of  a  variable  intrinsic 
impedance  for  specific  elemental  transmission  lines,  and  synchronism  is  maintained  by  increasing 
the  phase  velocity  along  the  diagonal  elemental  transmission  lines.  The  relationship  between  the 
FD  algorithm  and  TLM  model  is  established  through  the  equivalence  of  propagation  characteristics, 
the  most  fundamental  method  for  establishing  the  relationship  between  a  TLM  model  and  another 
numerical  method.  The  TLM  model  and  FD  algorithms  represent  identical  methods  for  the 
numerical  simulation  of  wave  propagation  if  the  latter  is  operated  at  the  upper  limit  of  its  stability 
range.  In  section  4.  the  propagation  characteristics  of  the  models  are  investigated.  For  the 
appropriate  selection  of  the  weighting  factor,  the  propagation  characteristics  become  approximately 
isotropic  (i.e..  the  directional  dependence  of  the  numerical  propagation  velocity  is  removed).  This 
allows  the  model  to  be  used  in  conjunction  with  the  velocity  error  correction  technique  described 
in  Reference  8.  In  section  5,  the  numerical  implementation  of  the  new  model  is  described.  The 
scattering  and  transfer  events  are  presented.  The  traditional  application  of  calculating  cutoff 
frequencies  in  rectangular  waveguide  is  used  to  validate  the  model.  Conclusions  and  a  discussion 
of  the  new  TLM  model  are  contained  in  section  6. 


2.  FINITE  DIFFERENCE  ALGORITHM 
Consider  the  following  semi-discretization  of  (2). 


EAx+M.v)  -  2E,(x.y)  -i-  £,(x-A/.y) 

£,(x.>+A/)-2£,(i,.v)  + £,(x.y-A/)  d^E, 

- - - (3) 

where  the  spatial  derivatives  are  replaced  with  second-order-accurate  central  difference  approxi¬ 
mations,  we  assume  cr  =  0.  and  the  right-hand  side  of  the  expression  is  evaluated  at  the  spatial 
location  (x.y).  The  stencil  for  this  spatial  discretization  is  shown  in  Figure  1.  We  assume  a  uniform 
grid  spacing  of  A/  in  the  x-  and  y-directions. 

FD  approximations  to  the  wave  equation  introduce  numerical  anisotropy  and  dispersion  (i.e.. 
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Figure  I.  Second-order-accuratc.  central  difference  Mcncil  used  to  approximate  the  «vpuiiai  derivatives  in  the  two-dimen* 

sional  wave  equation 


the  dependence  of  the  numerical  propagation  velocity  on  the  direction  of  propagation  and  fre¬ 
quency  content  of  the  signal).  To  reduce  the  numerical  aniswropy  present  m  the  semi-discretization 
(3),  Vichnevetsky  and  Bowles'"  proposed  the  weighted  combination  of  two  finite  difference 
approximations  to  the  spatial  derivatives  in  (2).  as  illustrated  in  Figure  2.  This  semi-discretization 
can  be  expressed  mathematically  as. 

(I-*)  E,{x-lLy) 

£,(x,y-t-d/)  -  2E.{x.y)  +  -A/)  1 

. .  A/-  '  I 

\E:(x+Al.y+SI)  -  2£.(x,v)  +  £.(x-A/.v-A/) 

I . . . UiAO: 

+  EAx+M.y-M)  -  2£,(x.>  )  -t-  £.(x-A/.v  +  A/)| 

(\2A/)-’  i 


Figure  2.  Weighted  combinalion  of  two  S-poini  stencils 
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where  A:  is  a  weighting  factor  restricted  between  zero  and  one  (and  again  the  right-hand  side  of 
(4)  is  evaluated  at  the  spatial  location  This  scheme  uses  the  same  grid  as  the  semi¬ 

discretization  (3)  and  has  the  same  memory  storage  requirements. 

With  the  appropriate  selection  of  k.  the  propagation  characteristics  of  (4)  become  isotropic  i.e.. 
the  propagation  velocity  becomes  independent  of  the  direction  of  propagation.  In  this  paper  we 
investigate  time-dependent  rather  than  time-harmonic  solutions  of  (2).  Therefore  approximation 
of  the  temporal  derivative  in  (4)  is  required.  Using  a  second-order-accurate  central  difference 
approximation,  (4)  becomes. 


{l-k) 


EAx+M.y)  -  lEAx.y)  +  £.(x-A/.y) 


£.(x.>'-t-4/)  -  2£.(x.v)  -t-  E.(x,y-M)] 

^  ^,2  - . 

^  ^|£,(x+A/.v+A/)  -_2EAx^)  ;t^EAx-M.y~ll) 
(^2A/)- 

^  ^(x  +  M.y-M)  -  2E,(x.y)  +  £,(x-AI.v  -(-A/) 
(\2MP 

E'/^lx.v)-2E'Ax.v)  +  E',  ^{x.v} 

- . . . - . 


where  At  denotes  the  time  step,  and  the  left-hand  side  of  (5)  is  evaluated  at  time  /.  (5)  represents 
an  explicitly  time-stepped  finite  difference  algorithm  for  the  solution  of  (2).  We  classify  this 
algorithm  as  an  explicitly-time-stepped,  second-order-accurate  in  time,  and  geometrically  weighted 
second-order-accurate  in  space,  FD  algorithm.  Trefethen"  has  investigated  this  algorithm  and 
determined  the  conditions  for  approximate  numerical  isotropy. 

The  dispersion  relation  for  a  numerical  method  yields  the  relationship  between  the  dispersed 
(or  numerical)  and  mathematically  exact  quantities.  We  use  the  notation  of  Vichnevetsky  and 
Bowles,'"  where  dispiersed  quantities  are  denoted  by  a(  *)  superscript  and  physical  (exact)  quantit¬ 
ies  are  otherwise  unscripted.  In  the  following  section,  the  dispersive  analysis  of  the  equivalent 
TLM  model  is  outlined.  It  is  necessary  to  distinguish  the  quantities  associated  with  the  elemental 
transmission  lines  of  the  model  from  both  the  numerical  and  physical  quantities.  We  use  an  (/) 
subscript  to  denote  elemental  transmission  line  quantities.  A  monochromatic  numerical  plane 
wave  propagating  through  the  numerical  mesh  at  an  angle  ih  to  the  x  axis  can  be  expressed  as. 


—  fi,  ^^'***  *■  10*  *  vsin*) 


(6) 


where  P*  represents  the  numerical  phase  constant.  Frequency  is  regarded  as  an  absolute  quantity 
defined  in  terms  of  numerical  or  exact  quantities. 


^  X*  '2Tr  X  2-n 


(7) 


where  c*  and  X*  are  the  numerical  propagation  velocity  and  wavelength,  respectively:  c.  X,  and 
p  are  the  exact  propagation  velocity,  wavelength  and  phase  constant,  respectively  (c  =  (cp.)"''^). 
Substitution  of  (6)  into  (5)  yields  the  disfiersion  relation  for  the  finite  difference  algorithm. 


f  -  -P* 

(sin*  — 


A/(cos  <{>  +  sin  <|>)  .  ,p*A/(cos<l)  -  sin<l>) 

- - +  sin-  — - - - 


+  (I-k) 


.  ,p*A/cos<l>  .  ,p*A/sin«f> 
sin* - - - +  sin*  - - 


A/-  .  ,  loA/ 

,T-^  sin* 
c*Ar-  2 


(8) 


Expression  (8)  describes  the  fundamental  manner  in  which  plane  waves  propagate  through  an 
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infinite  FD  grid.  Given  a  spatial  and  tempoial  discretization  (A/  and  Ar  respectively),  frequency 
of  excitation  (oi),  direction  of  propagation  (<j>).  and  the  weighting  factor  (^),  the  numerical  phase 
constant  (3*)  can  be  obtained  from  (8).  This  value  can  be  compared  to  the  exact  physical  phase 
constant  to  determine  the  amount  of  dispersion  introduced  by  the  algorithm.  Therefore.  (8)  is  a 
fundamental  representation  of  the  fidelity  of  the  algorithm  as  a  method  for  the  simulation  of  wave 
propagation. 

TTie  stability  criterion  for  this  algorithm  (obtained  using  the  Von  Neumann  method,  discussed 
in  Reference  12)  is  given  by. 


Ars 


I  ^ 
,2^  ^ 


(9) 


3.  TRANSMISSION-LINE  MATRIX  MODEL 


3.1.  Synthesis 

We  now  synthesize  a  TLM  model  equivalent  to  the  FD  algorithm  presented  in  the  previous 
section.  The  FD  algorithm  is  constructed  from  the  weighted  combination  of  two  second-order- 
accurate  central  difference  stencils,  one  oriented  as  usual  (arms  of  the  stencil  located  along  the  x 
and  y  axis),  the  other  rotated  by  45"  with  its  arms  extended  by  a  factor  of  (2)'  *.  It  has  been 
demonstrated  that  the  original  TLM  model’  and  the  FD  algorithm  (3)  (with  temporal  derivatives 
approximated  by  a  second-order-accurate  central  difference  approximation )  are  e  ,uivalent.*  There¬ 
fore,  the  new  TLM  model  should  consist  of  the  weighted  combination  of  two  original  models. 
One  oriented  as  usual  (with  elemental  transmission  lines  oriented  along  the  x  and  ,v  axis),  the 
other  rotated  by  45°  with  its  arms  extended  by  a  factor  of  (2)*'^.  The  basic  geometry  of  the  model 
is  shown  in  Figure  3.  The  new  model  is  realized  as  a  shunt  connection  of  transmission  lines  (as  in 
Reference  1).  A  mesh  of  nodes  is  provided  in  Figure  4.  Note  that  a  direct  electrical  connection 
between  diagonal  and  axial  transmission  lines  exists  only  at  the  centres  of  nodes,  located  at  even 
multiples  of  A/  in  both  the  x-  and  y-directions  (denoted  by  the  black  dots  in  the  figure).  To 


Figure  3.  Basic  geometry  of  the  new  TLM  model,  created  from  the  combination  of  two  original  models 
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Figure  4.  Mesh  of  new  nodes.  Electrical  connections  in  the  mesh  are  denoted  by  black  dots.  These  spatial  locations  are 

the  centres  of  nodes 


complete  the  model,  the  electrical  characteristics  of  the  elemental  transmission  lines  must  be 
determined. 

The  electrical  circuit  analogue  of  a  weighting  factor  is  a  variable  impedance.  To  implement  the 
ability  to  weight  the  two  interconnected  original  models,  the  diagonal  and  axial  elemental  trans¬ 
mission  lines  are  permitted  to  have  different  characteristic  impedances.  The  intrinsic  impedance 
of  the  axial  elemental  transmission  lines  (i.e.,  associated  with  the  original  model  with  elemental 
transmission  lines  along  the  x  and  y  axis)  is  Z,,  and  the  intrinsic  impedance  of  the  diagonal 
transmission  lines  (i.e..  associated  with  the  original  model  rotated  by  45°)  is  mZ,.  where  m  is  the 
impedance  weighting  factor  (0  «  m  <  infinity). 

In  the  evaluation  of  the  FD  algorithm  (5),  communication  of  information  between  spatial 
locations  in  the  axial  direction  takes  place  at  the  same  speed  as  communication  of  information 
between  spatial  locations  in  the  diagonal  direction.  Therefore,  propagation  along  diagonal  elemen¬ 
tal  transmission  lines  should  be  (2)*'^  times  faster  than  propagation  along  the  axial  elemental 
transmission  lines,  or 


V';=y2v;  (10) 

where  v"  refers  to  the  propagation  velocity  along  the  nth  elemental  transmission-line,  j  =  5-8, 
and  i  =  1-4.  A  beneficial  consequence  of  (10)  is  that  the  synchronism  of  voltage  pulses  is 
preserved  in  the  new  model.  The  electrical  and  geometrical  description  of  the  new  model  is 
complete. 


3.2.  Propagation  analysis 

The  topology  of  the  model  is  provided  in  Figure  5.  To  model  a  medium  of  arbitrary  permntivity, 
an  open  circuit  stub  is  added  to  the  centre  of  a  TLM  node.'*  The  new  model  is  the  weighted 
combination  of  two  original  shunt  nodes.  Therefore,  to  maintain  consistency,  two  open  circuit 
stubs  are  added.  One  of  length  A//2  and  admittance  YolZi  (associated  with  the  shunt  node  with 
elemental  transmission  lines  along  the  x  and  y  axis),  the  other  length  and  admittance 

YglmZi  (associated  with  the  rotated  shunt  node). 

The  propagation  analysis  of  the  model  proceeds  in  the  same  manner  as  performed  in  References 
6,  8  and  13.  Superposition  and  transmission-line  theory  yield  the  characteristic  equation  which 
describes  the  behaviour  of  voltages  on  the  model, 

2(m+ 1 )  [2  cos  ^Al  -  F()  sin^  v,  +  Vj 


(11) 
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Figure  5.  Complete  electrical  and  geometnc  description  of  the  new  model  with  permittivity  stubs 


A  monochromatic  numerical  plane  voltage  wave  propagating  through  ihe  numerical  mesh  at  an 
angle  4)  to  the  x  axis  can  be  expressed  as. 


V 


V„e 


fUt/  ^  |(i“  (  V  ♦  X  Mf)  tl>) 


(12) 


where  the  parameters  in  the  exponential  term  of  (12)  are  as  defined  in  section  2.  Substitution  of 
(12)  into  (11)  yields  the  dispersion  relation  for  the  TLM  model. 


.  ,  p*A/(cos<t)  +  sin<i))  .  ,  3*A/(cos<i>  -  sin  di) 
sin-  - — — r - +  sin-  - ^ - 


+  m  isin 


,P*A/cos<t>  .  ,p*A/sin(l>]  w-i  .  , .  ,  .  ,3/A/ 

^2  — - - +  - - —  i  =  (4+  5,„. 


(13) 


Expression  ( 13)  describes  the  fundamental  manner  in  which  plane  waves  propagate  through  an 
infinite  TLM  mesh.  Given  a  spatial  discretization  (A/),  frequency  of  excitation  (described  through 
B/),  direction  of  propagation  (4)),  and  the  electrical  properties  of  the  model  (m  and  y',,)-  the 
numerical  phase  constant  (3*)  can  be  obtained  from  (13).  This  value  can  be  compared  to  the 
exact  physical  phase  constant  to  determine  the  amount  of  dispersion  introduced  by  the  model. 
Therefore,  (13)  is  a  fundamental  representation  of  the  fidelity  of  the  model  as  a  method  for  the 
simulation  of  wave  propagation. 


3.3.  Equivalence  of  the  TLM  model  and  FD  algorithm 

We  now  establish  the  equivalence  of  the  TLM  and  FD  methods  and  demonstrate  that  both  can 
represent  identical  models  for  wave  propagation.  This  is  accomplished  by  determining  the  con¬ 
ditions  for  which  (8)  and  (13)  are  equivalent. 

The  term  3/^/  in  the  right-hand  side  of  (13)  can  be  re-expressed  as  u>Af  by  noting  the  following 
relationships, 


2'?r 

(14a) 

2irv; 

X,= - ? 

(1> 

(14b) 

M 

‘’'=Ar 

(14c) 

(14b)  is  a  direct  extension  of  (7),  (i.e..  frequency  is  considered  as  an  absolute  quantity  and  can 
be  defined  in  terms  of  exact,  numerical,  or  elemental  transmission-line  quantities).  If  we  divide 
the  FD  dispersion  relation  (8)  by  kl2  we  obtain. 
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.  ,  p*A/(cos4)  +  sin<|>)  .  ,  3'A/(cos<i)  -  siinji) 

,n-  - - - +  sin- - - — - - 


sin 


2(1-^) 


f  .  ,  S*A/  cos  (J) 

— 


sin‘ 


,3»A/sin<i.]  21F 


.  ijjAt 


(15) 


(15)  and  (13)  have  fundamentally  the  same  form.  Equating  coefficients  of  the  left-hand  side  of 
(8)  and  (13)  yields  the  relationship  between  the  TLM  and  FD  weighting  factors, 


m  = 


2(1-^) 


(16a) 


or 


m  +  2 


(16b) 


Equating  the  coefficients  on  the  right-hand  side  of  ( 15)  and  ( 13)  yields. 


c  = 


_ 2 _ A/ 


(I7a) 


substitution  of  (I6a)  into  (17a)  yields. 


r=  2  A/ 

v(2-A:)(4+K„)-^' 


(17b) 


or  if  we  desire  c  in  terms  of  TLM  model  parameters  alone,  substitution  of  ( I6a)  and  ( 14c)  into 
( 17b)  yields. 


j  2(^2) 

^  V(m+l)(4-(-K.,)‘’' 


(17c) 


If  the  FD  algorithm  is  operated  such  that  (17b)  is  satisfied,  the  dispersion  relations  for  both  are 
identical,  and  therefore  the  two  methods  fundamentally  represent  identical  methods  for  the 
simulation  of  wave  propagation. 

It  is  interesting  to  note  that  for  the  condition  T,,  =  0  (a  ‘free  space’  TLM  model),  the  condition 
(17b)  corresponds  to  the  upper  limit  of  the  FD  stability  range.  As  was  found  for  the  original  node 
and  the  Yee  algorithm,*  the  TLM  model  and  FD  algorithm  are  identical  when  the  latter  is  operated 
at  the  upper  limit  of  its  stability  range.  This  was  not  the  case  for  the  hexagonal  TLM  and  FD 
methods.** 

If  we  return  to  the  context  of  modelling  electromagnetic  phenomena,  we  can  establish  the 
relationship  between  the  admittance  of  the  open  circuit  stub  and  the  material  properties  of  the 
medium  modelled  by  the  entire  model.  The  physical  propagation  velocity  is  defined  as. 


1 

c  = . . . 


(18) 


where  e,  and  p,,  are  the  relative  permittivity  permeability,  respectively  and  Co  and  p,)  are  the  free 
space  permittivity  and  permeability,  respectively.  Relating  this  to  the  propagation  velocity  in  the 
TLM  model  (given  by  (17c)),  we  obtain. 


1  j  2{m+2)  M 

vw:;:;^V(m+i)(4+y„)At 


(19) 


If  we  consider  the  case  Y„  =  0  to  represent  free  space,  i.e.,  e,  =  p,  =  1,  (19)  becomes. 
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1  ^  lm+2)  AP 

Cofio  2(m-t-l) 

The  relative  permittivity  and  relative  permeability  are  related  to  the  stub  admittance  T,,  by, 

e.iL,  =  (l  +  y)  (21) 


4.  PROPAGATION  CHARACTERISTICS 

Numerical  models  for  wave  propagation  represent  a  discretized  medium  that  is  both  dispersive 
and  anisotropic,  i.e.,  the  propagation  velocity  of  waves  in  the  numerical  mesh  depends  on  both 
the  frequency  content  of  a  signal  and  the  direction  of  propagation.  This  undesired  effect  is  referred 
to  as  velocity  error  and  is  determined  from  the  dispersion  relation  for  the  particular  model.  The 
ratio  c*/c  (the  ratio  of  the  numerical  propagation  velocity  to  the  physical  propagation  velocity), 
can  be  used  as  a  quantitative  measure  of  velocity  error.  The  TLM  dispersion  rela  ion  ( 13)  can  be 
rewritten  as. 


,  -il  ,  A/ 

sin*  n-fcGsd  +  sin  4>)  sin-  •ir(cos  d)  —  sin  d))  ; 

A  \ 


+  m 
m+l 


^  A/  ,  .  ,  ^  A/1 

sin-  TT  cos  4>  r;  sin-  it  cos  d>  1 
A  X 


(4+ y,,)  sin- IT 


2(m+2)  A/ 


(w-t-lKTo-^-A)  X 


(22) 


Given  the  free  space  discretization  ratio  (A//X).  direction  of  propagation  (6).  and  the  electrical 
properties  of  the  model  (m  and  Y'»),  (22)  can  be  searched  to  determine  the  dispersed  discretization 
ratio  (A//X*).  Given  A//X  and  A//X*.  the  ratio  c*/c  can  be  determined  from. 


c*  _  A//X 
c  ■■  A//X* 


(23) 


In  Figure  6(a),  (b),  (c),  (d),  (e),  and  (f),  c*/c  is  provided  as  a  function  of  d)  for  m  =  900,  6,  4. 
3,  2  and  0  01,  respectively.  For  each  case,  contours  for  A//X  =  OTO.  0-20.  0-30.  and  0-35  are 
provided  (To  =  0  for  all  cases).  Note  that  in  light  of  the  equivalence  established  in  section  3.3, 
Figure  6(a),  (b),  (c),  (d),  (e)  and  (f)  are  applicable  to  the  FD  algorithm  provided  the  FD 
algorithm  is  operated  at  its  upper  limit  of  stability,  c,  =  M-r  =  1  and  k  =  0  002217.  0-25,  0-333. 
0-4,  0-5,  and  0-995,  respectively. 

In  the  limit  as  m  approaches  infinity,  the  new  model  is  equivalent  to  a  mesh  of  original  nodes’ 
with  elemental  transmission  lines  oriented  along  the  x  and  y  axis.  In  Figure  7.  c*/c  is  provided  as 
a  function  of  d>  for  the  original  model'  (for  F,,  =  0).  As  expected  the  contours  of  Figure  7  and 
Figure  6(a)  are  indistinguishable.  In  the  limit  as  m  approaches  zero,  the  new  model  is  equivalent 
to  the  original  model  rotated  by  45°  and  mesh  spacing  extended  by  a  factor  of  (2)'  -.  In  Figure 
8,  c*/c  is  provided  as  a  function  of  tb  for  the  original  model  rotated  by  45°  and  mesh  spacing 
extended  by  a  factor  (2)''^  (for  F,,  =  0).  As  expected  the  contours  of  Figure  8  and  Figure  6(f) 
are  indistinguishable. 

For  moderate  values  of  m,  directions  for  propagation  with  no  dispersion  do  not  exist  with  the 
new  model.  From  the  results  of  Figure  7.  we  note  that  no  numerical  dispersion  exists  for  waves 
which  propagate  diagonally  through  the  mesh  (<}>  =  45°  +  n90°,  n  =  0,  1.  2.  3).  Numerical 
dispersion  is  maximum  for  axial  propagation  (ifc  =  n90°,  n=0.  1,  2,  3).  For  the  rotated  original 
model,  the  complementary  situation  is  present.  No  numerical  dispersion  exists  for  (di  =  «90°.  n 
=  0,  1,  2,  3),  and  numerical  dispersion  is  maximum  for  (4)  =  45°  +  n90°.  n  =  0,  1.  2.  3).  From 
Figure  6(b)-(e)  we  note  that  the  new  model  blends  the  propagation  characteristics  of  the  original 
and  rotated  original  models.  Therefore,  propagation  along  the  directions  for  maximum  numerical 
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Figure  6.  Contours  of  c*/c  for  the  new  TLM  modcJ  (Vo  =  0)  for  (a)  m  9(X),  (b)  m  =  6.  (c)  m  -  4.  (d)  m  »  3.  (c) 

m  -  2  and  (f)  m  -  0  01 


dispersion  is  improved,  but  directions  for  perfect  propagation  are  eliminated.  Therefore,  in  this 
context  the  propagation  characteristics  of  the  original  model  are  superior  to  those  of  the  new 
model. 

However,  from  Figure  6  it  can  be  noted  that  for  the  appropriate  selection  of  the  weighting 
factor,  the  new  model  can  possess  propagation  characteristics  with  approximate  isotropy.  The 
appropriate  conditions  have  been  investigated  in  the  context  of  the  equivalent  FD  algorithm.'”  " 
TTie  appropriate  weighting  factor  for  the  semi-discretization  (4)  is  it  =  ()-5  (see  Reference  lO). 
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Figure  6.  Continued 


For  the  full  discretization  (5),  Trefethcn"  has  determined  a  weighting  factor  of  it  =  1/3  provides 
isotropy  to  order  (A/)*  (note  that  the  difference  in  k  for  the  semi-discretization  and  full  discretiz¬ 
ation  is  a  result  of  effect  of  temporal  discretization  in  the  later).  Therefore,  isotropy  to  order 
(A/)*  should  be  obtained  from  the  new  TLM  model  for  m  =  4  0  (using  (16a)  to  convert  k  to  m). 
as  shown  in  Figure  6(c). 

Obtaining  approximate  numerical  isotropy  is  equivalent  to  reducing  the  dependence  of  the 
propagation  velocity  on  the  direction  of  propagation.  Consider  the  simulation  of  a  homogeneous 
problem  that  employs  a  regular  mesh.  If  the  numerical  propagation  velocity  is  independent  of  the 
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Figure  7.  Contours  of  c“‘c  for  the  original  TLM  model  ‘  Contours  are  mdistinguishabic  from  ihiKc  of  Figure  Wa) 
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Figure  8.  Contours  of  c*  c  for  the  original  TLM  model  rotated  by  45®  and  mesh  spacing  extended  by  a  factor  (2)‘  *. 
Contours  are  indistinguishable  from  those  of  Figure  6(f) 
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direction  of  propagation,  the  amount  of  velocity  error  at  a  given  frequency  can  be  accurately 
estimated  from  the  dispersion  relation.  Therefore,  after'the  simulation  is  complete,  the  velocity 
error  can  be  corrected  at  each  output  frequency.  This  method  was  used  in  Reference  8  (with  the 
hexagonal  two-dimensional  TLM  model)  to  correct  for  the  cutoff  frequencies  in  a  rectangular 
waveguide  with  no  a  priori  assumption  regarding  the  directional  dependence  of  a  particular  mode. 
In  this  context  the  propagation  characteristics  of  the  new  model  can  be  considered  as  superior  to 
those  of  the  original  model.' 

In  Figure  9,  the  propagation  characteristics  of  the  new  model  and  the  hexagonal  TLM  model* 
are  compared.  The  ratio  c*/c  is  provided  as  a  function  of  the  physical  discretization  ratio  (A//X) 
for  propagation  directions  ih  =  0°,  22-5°.  39°,  and  45°  for  (a)  the  new  mode!  with  m  =  3-0,  (b) 
the  new  model  with  m  =  4  0  and  (c)  the  hexagonal  TLM  model.  The  results  contained  in  the 
figure  indicate  that  the  hexagonal  model  is  superior  to  the  new  model  in  terms  of  both  the  cutoff 
frequency  of  the  model,  and  the  degree  of  approximate  isotropy.  Therefore,  in  the  context  of 
isotropic  models  for  the  simulation  of  wave  propagation,  the  hexagonal  model  is  preferred. 

An  advantage  of  the  new  node  is  that  it  is  realized  on  a  regular  grid  with  equal  spacing  in  the 
X-  and  y-directions  (A/)  The  hexagonal  model  is  also  realized  using  a  mesh  with  equal  inter-nodal 
spacing.  However,  owing  to  the  nature  of  the  hexagonal  grid,  the  spacing  in  the  x-  and  y-directions 
is  unequal,  A/  and  (3)‘'^A//2,  respectively.  This  creates  a  disadvantage  for  the  hexagonal  model 
in  the  modelling  of  structures  with  regular  geometric  features. 
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Figure  9.  c*/c  versus  A//\  (the  exact  (physical)  discretization  ratio)  for  vanous  of  <f>  for  (a)  the  new  model  with  m  =>  3. 
(b)  the  new  model  with  m  ^  4  and  (c)  the  hexagonal  TLM  model 

5.  NUMERICAL  IMPLEMENTATION  AND  RESULTS 


5.1.  Scattering  matrix 

In  the  previous  sections  we  have  examined  some  of  the  theoretical  aspects  of  the  new  TLM 
model.  We  now  describe  the  numerical  implementation  of  the  model  in  terms  of  the  traditional 
scattering  and  transfer  events.  TLM  algorithms  operate  by  simulating  the  progression  of  voltage 
pulses  as  they  are  scattered  throughout  the  mesh  of  transmission  lines.  Applying  the  appropriate 
initial  conditions  and  reflection  coefficients  (to  model  boundary  conditions)  the  transmission-line 
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simulation  becomes  analogous  to  an  electromagnetic  field  problem.  The  implementation  of  the 
new  model  follows  the  same  procedure  as  all  other  models,  i.e..  scattering  of  incident  impulses 
at  the  junction  of  nodes  and  transfer  of  the  reflected  pulses  to  adjacent  nodes.  The  algorithm  can 
be  expressed  formally  as. 


,V'  =  S,V‘  (24) 

and 

,.,V‘  =  C,V'  (25) 

where  and  are  the  vectors  of  the  incident  and  reflected  voltage  pulses  at  all  nodes  at  time- 
step  k.  S  is  the  global  scattering  matrix  describing  the  interaction  of  pulses  at  all  nodes  in  the 
mesh,  and  C  is  the  connection  matrix  describing  how  nodes  are  connected  (and  includes  the 
boundary  conditions  for  the  particular  problem).  These  two  equations  include  all  information 
required  to  perform  the  simulation. 

The  nodal  scattering  matrix  can  be  assembled  by  examining  the  reflection  and  transmission 
coefficients  of  a  voltage  pulse  on  each  of  the  ten  elemental  transmission  lines  of  the  model.  A 
voltage  pulse  on  the  /th  elemental  transmission  line  sees'  a  reflection  coefficient  of. 


,--2'  - 

z,.  +  z\ 


(26) 


where  Z^.  is  the  parallel  combination  of  all  but  the  ith  elemental  transmission  line  and  Z)  is  the 
intrinsic  impedance  of  the  ith  elemental  transmission  line.  The  intrinsic  impedance  of  the  elemental 
transmission  lines  (from  section  3.1  and  shown  in  Figure  5).  is. 


z;- 


Z/  for  I  =  1-4 

mZ,  for  /  =  5-H 

Z,/Y„  for  I  =  d 

mZilYo  for  /  =  10 


The  associated  transmission  coefficient  is. 


r=  1  +  r 

From  (26)-(28),  the  nodal  scattering  matrix  can  be  assembled  as. 
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=  vj{/-l.y) 

‘'3  (ij)  =  *'{  (».y>  1 ) 
vkU.j)  =  v':ii+l.i) 
y'sU^i)  =  vni-i./-i) 

K{iJ)  =  vj(/-l./+l) 
FtCi./)  =  v?(i+l,/+l) 

=  V6(i+i.;-i) 

V‘<>  UJ)  =  (i,/) 

‘'io  (<7)  =  Vio  ('.» 


(30) 


where  (i,/)  denotes  the  discrete  location  of  a  node  in  the  mesh. 

If  the  TLM  method  is  considered  as  a  differential-equation-based  numerical  method  for  solving 
(2),  (29)  and  (30)  represent  the  approximate  model  for  wave  propagation  (in  the  same  way  the 
FD  method  is  considered  as  a  differential-equation-based  numerical  method  for  solving  (2)  and 
(S)  represents  the  approximate  model  for  wave  propagation).  The  solution  of  a  specific  problem 
requires  the  application  of  initial  and  boundary  conditions.  The  treatment  of  boundary  conditions 
is  an  important  subject  for  the  practical  application  of  the  method.  In  this  paper  we  are  primarily 
concerned  with  the  development  of  the  new  TLM  model  as  an  approximate  model  fot  ‘'’ave 
propagation  and  establishing  the  equivalence  with  the  FD  algorithm.  Therefore  we  do  not  treat 
the  subject  of  boundary  conditions  in  detail.  The  traditional  methods  of  specifying  reflection 
coefficients  at  locations  half-way  between  the  centres  of  nodes  (in  both  the  axial  and  diagonal 
direction)  should  be  applicable. Potential  users  should  be  cautioned  that  the  intrinsic  impedance 
of  the  elemental  transmission  lines  is  not  always  the  same  for  this  model  and  care  should  be  taken 
in  the  evaluation  of  the  appropriate  reflection  coefficients  for  a  specific  boundary  condition.  The 
method  described  by  Chen  ei  al.'*  of  enforcing  boundary  conditions  at  the  centre  of  nodes  should 
also  be  applicable  to  the  new  model. 


5.2.  Calculation  of  cutoff  frequencies 

To  validate  the  new  TLM  model,  we  investigate  the  traditional  TLM  application  of  the  calcu¬ 
lation  of  cutoff  frequencies  of  various  modes  in  a  waveguide.^  ^  The  cross-section  of  the  partially 
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Table  I.  Cutoff  frequencies  (in  GHz)  for  rectangular  waveguide 

Mode 

New  TLM  model 

Analytical'' 

“o  difference 

(1.1) 

1-244 

1-2490 

0-4 

(1.2) 

1-784 

1-8015 

0-9 

(2.1) 

2  126 

2-1345 

0-4 

(1.3) 

2-429 

2-4605 

1-3 

(2.2) 

2-481 

2-4983 

0-7 

Table  11. 

Cutoff  frequency  (in  GHz)  of  the  dominant  mode  in  par- 

tially  filled  waveguide 

«i 

New  TLM  model 

Finite  element"’ 

“o  difference 

2-5 

1-054 

1  063 

()-8 

50 

0-846 

0-852 

0-7 

100 

0614 

0-623 

1-4 

filled  rectangular  waveguide  and  the  physical  dimensions  are  provided  in  Figure  10.  The  walls  of 
the  guide  are  considered  to  be  perfectly  conducting.  To  realize  this  boundary  condition,  reflection 
coefficients  of  magnitude  -10  are  placed  at  locations  half-way  between  nodes.  A  mesh  spacing 
of  Al  =  0  01  metres  is  selected,  resulting  in  a  total  TLM  mesh  with  20  nodes  in  the  jr-direction 
and  15  nodes  in  the  y-direction.  Calculations  were  performed  such  that  the  true  physical  cutoff 
frequencies  are  obtained  directly  from  the  simulation.  Normalization  for  a  non-free  space  medium 
is  not  required. 

The  TLM  simulation  yields  the  cutoff  frequencies  of  TM  modes.  Table  1  contains  the  cutoff 
frequencies  for  the  first  five  modes  for  c,  =  10  {¥„  =  0).  A  total  of  1000  iterations  (i.e..  lOOOAr, 
where  At  can  be  obtained  from  (17c))  and  a  weighting  factor  of  m  =  6  was  used.  The  TLM  results 
are  compared  to  analytical  results.'*  The  percentage  difference  is  provided  in  the  table.  Reasonable 
accuracy  is  obtained.  In  Table  11.  the  cutoff  frequency  of  the  dominant  mode  is  provided  for  c, 
=  2  5,  5  and  10.  As  a  comparison  results  generated  by  a  finite  element  code  are  provided.'* 
Again,  reasonable  agreement  is  obtained. 


6.  CONCLUSIONS 

In  this  paper  we  have  presented  a  new  TLM  model  for  the  simulation  of  the  two-dimensional 
wave  equation.  The  TLM  model  was  synthesized  directly  from  an  FD  algorithm"’ "  as  a  shunt 
connection  of  two-wire  transmission  lines.  The  new  model  is  a  spatially  weighted  connection  of 
two  original  models.'  One  oriented  as  usual,  the  other  rotated  by  45°.  The  weighting  is 
accomplished  through  the  use  of  a  variable  intrinsic  impedance  for  specific  elemental  transmission 
lines.  Synchronism  is  maintairicd  by  increasing  the  propagation  velocity  along  diagonal  elemental 
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transmission  lines.  In  general,  the  synthesis  of  other  TLM  models  from  FD  algorithms  is  possible 
using  the  same  basic  steps. 

The  relationship  between  the  FD  algorithm  and  TLM  model  is  established  through  the  equival¬ 
ence  of  propagation  characteristics.  We  feel  this  is  the  most  fundamental  method  for  establishing 
the  relationship  between  a  TLM  model  and  another  numerical  method.  It  is  possible  to  demonstrate 
that  the  TLM  model  satisfies  the  FD  algorithm  (by  examining  the  scattering  and  transfer  of  voltage 
pulses  at  a  node  and  its  neighbours).  In  Reference  17.  the  original  TLM  model  was  shown  to 
satisfy  the  two-dimensional  Yee  algorithm.  To  demonstrate  the  equivalence,  the  definitions  for 
magnetic  field  quantities  in  terms  of  voltage  pulses  were  altered.  Rather  than  define  the  magnetic 
field  components  at  the  centres  of  nodes  at  each  iteration,  the  magnetic  fields  were  defined  at  the 
intersection  of  nodes  at  half  iterations.  Chen  ei  al.'*  have  recently  reported  a  complete  algorithmic 
equivalence  of  two-  and  three-dimensional  TLM  models  and  FD  algorithms.  We  regard  the 
calculation  of  field  values  from  voltage  pulses  to  be  a  post-processing  task  associated  with  the  TLM 
method.  The  basic  algorithm  (scattering  and  transfer  of  voltage  pulses)  operates  independently  of 
these  definitions.  One  of  the  often-quoted  advantages  of  the  TLM  approach  is  the  ability  to  define 
field  compionents  at  various  spatial  and  temporal  locations  (as  long  as  a  certain  consistency  is 
maintained).  Therefore,  we  feel  that  establishing  an  equivalence  between  a  TLM  model  and 
another  numerical  method  without  specific  definitions  for  field  quantities  in  the  TLM  model, 
is  the  most  fundamental  and  rigorous.  We  accomplish  this  by  demonstrating  the  propagation 
characteristics  of  the  TLM  model  and  FD  algorithm  are  identical  if  the  latter  is  operated  at  the 
upper  limit  of  its  stability  range. 

The  propagation  characteristics  of  the  new  model  have  been  examined.  For  moderate  values  of 
m.  directions  for  propagation  with  no  dispersion  do  not  exist  for  the  new  model.  Therefore,  in 
this  context  the  propagation  characteristics  of  the  original  model  are  superior  to  those  of  the  new 
model.  However,  an  advantage  of  the  new  model  is  that  for  appropriate  values  of  m,  approximate 
numerical  isotropy  is  obtained.  This  allows  the  model  to  be  combined  with  an  error-correction 
method  to  remove  the  contribution  of  velocity  error  from  the  results.*^  Comparison  of  the  character¬ 
istics  of  the  new  model  to  those  of  the  hexagonal  TLM  model'*  indicate  the  hexagonal  model  is 
preferred  (in  terms  of  both  the  amount  of  approximate  numerical  isotropy  and  cutoff  frequency). 
This  conclusion  is  supported  by  analogous  finite  element  (FE)  studies.  Consider  the  relationship 
of  the  various  FD  algorithms  and  TLM  models  ( References  6  and  8.  and  section  3  of  this  paper), 
and  the  relationship  of  the  FD  and  FE  methods.*'' Based  on  these  relationships,  the  original 
TLM  model*  is  analogous  to  using  square  quadrilateral  finite  elements  of  sides  ,>/;  the  hexagonal 
TLM  model**  is  analogous  to  using  equilateral  triangular  finite  elements,  each  triangle  having  sides 
of  A/  and  angles  of  60°;  and  the  new  TLM  model  is  in  some  way  analogous  to  using  right  triangular 
finite  elements,  each  triangle  having  sides  of  A/.  A/.  (2)'  -A/  and  angles  90°.  45°,  45°.  Mullen  and 
Belytschko  have  determined  that  modelling  with  equilateral  triangles  ( analogous  to  the  hexagonal 
TLM  model)  is  the  optimum  triangular  discretization  if  isotropy  is  desired.-*  This  supports  the 
analysis  performed  in  this  paper. 

Finally,  the  scattering  and  transfer  events  for  the  new  model  were  presented  and  were  applied 
to  the  analysis  of  a  rectangular  waveguide  partially  filled  with  a  dielectric.  The  cutoff  frequencies 
calculated  using  the  new  TLM  model  agreed  well  with  both  analytic  and  numerical  finite  element 
results. 

In  1976.  Johns  presented  an  interesting  papier  in  which  the  original  TLM  model  is  described  as 
a  discrete  form  of  Huygens'  Principle.--  Hoefer  has  continued  this  view  and  has  provided  a  brief 
historical  review  and  description  of  the  discretization  process.^  The  hexagonal  model  can  be 
considered  as  a  logical  extension  of  the  original  model.  The  improvement  in  numerical  isotropy 
over  the  original  model  is  intuitively  obvious.  The  model  presented  in  this  paper  could  also  be 
described  as  a  discrete  form  of  Huygens'  Principle.  However,  the  model  would  have  been 
developed  with  a  specific  value  for  the  weighting  factor  (presumably  such  that  energy  would  be 
scattered  isotropically).  While  selecting  a  variable  weighting  factor  may  be  of  more  theoretical 
interest  than  practical  value,  the  motivation  for  allowing  this  flexibility  may  not  be  obvious  from 
the  perspective  of  a  discrete  form  of  Huygens'  Principle, 

The  original  model,'  the  hexagonal  model*  and  the  new  model  presented  in  this  paper  are 
equivalent  to  FD  algorithms  that  approximate  spatial  derivatives  with  second-order-accurate  central 
d'Terence  formulas.  The  difference  between  the  various  models  is  the  geometric  configuration 
and  weighting  of  the  difference  approximations.  Future  work  will  investigate  the  synthesis  of  a 
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TLM  model  equivalent  to  an  FD  algorithm  that  approximates  the  spatial  derivatives  in  the  wave 
equation  with  fourth-order-accurate  central  difference  formulas. 
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SUMMARY 

The  state  and  output  equations  of  the  overall  networks  are  derived  from  the  state  and  output  equations  of 
individual  multiports  and  knowledge  of  the  interconnections  between  them.  A  generalized  lumped-distributed 
L/D  multiport  is  described  by  its  associated  state,  output  and  non-linear  equations  in  the  time  domain.  Any 
network  can  be  considered  as  composed  of  a  set  of  multiports  and  independent  sources.  These  equations 
have  been  incorporated  into  a  computer-aided  procedure  for  the  analysis  of  L/D  networks.  The  procedure 
can  be  used  for  the  simulation  of  any  non-linear  microwave  circuit  and  offers  the  facility  of  developing  a 
multiport  equivalent  circuit  for  any  linear  or  non-linear  device  or  subcircuit.  Several  examples  are  successfully 
analysed  using  the  developed  general  program. 

1.  INTRODUCTION 

The  analysis  of  non-linear  dynamic  networks  by  using  state-space  approach  has  been  established 
since  the  1960s  and  well  documented  in  many  reference  books.*  *  Computer-aided  state-space 
analysis  of  lumped  and  lumped/distributed  networks  has  been  developed.'  The  capacitor 
voltages  (or  charges)  and  the  inductor  currents  (or  fluxes)  are  usually  chosen  as  the  lumped  state 
variables.  The  reflected  voltages  at  the  transmission  lines  can  be  chosen  as  the  distributed  state 
variables.  In  all  these  cases  the  state  and  output  equations  are  established  from  the  circuit  element 
values  and  the  topology  of  the  whole  network. 

It  is  highly  desirable  and  convenient  for  circuit  designers  to  consider  the  non-linear  network 
composed  of  subcircuits.  These  subcircuits  are  represented  by  functional  blocks  described  by  a  set 
of  equations.  In  this  case  the  formulation  of  the  whole  network  equations  starts  from  the  top  level 
of  the  subcircuits  (multiports).  The  graph  of  the  network  is  only  describing  the  interconnection 
of  all  network  multiports. 

Multiport  representation  is  common  for  linear  networks  in  the  frequency  domain  where  any  of 
the  usual  multiport  parameters  (x,  y,  h,  g)  or  the  scattering  parameters  can  be  used.  When  any 
of  these  parameters  are  known,  the  topology  and  element  values  of  the  multiport  are  no  longer 
required.  No  such  treatment  has  so  far  been  available  for  non-linear  networks.  Non-linear  networks 
are  usually  solved  in  the  time  domain  either  by  direct  integration  of  the  network  equations,'  by 
using  associated  discrete  circuit  modelling  (Spice)  or  by  the  harmonic  balance  technique.^  In  the 
harmonic  balance  method  the  non-linear  subnetwork  is  still  solved  in  the  time  domain.  It  is  then 
a  great  advantage  to  develop  a  method  of  characterizing  non-linear  networks  from  their  terminal 
behaviour  and  treat  them  as  multiports. 

In  this  work  multiports  can  represent  networks  with  I  imped,  distributed  and  non-linear  elements. 
Each  multiport  is  represented  by  non-linear  state  and  output  equations  and  the  overall  network 
is  composed  of  a  number  of  individual  multiports  connected  in  an  arbitrary  fashion.  The  state  and 
output  equations  of  the  overall  network  are  derived  and  solved  in  the  time  domain.  Thus  the 
method  enables  the  hierarchical  development  of  non-linear  networks.  At  the  lowest  level  of 
hierarchy  the  multiport  equivalent  is  developed  from  individual  circuit  elements  (linear  and/or 
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non-linear)  using  a  tabular  approach.  At  the  higher  levels  only  the  multiport  equivalent  is  required. 
Any  number  of  hierarchical  levels  can  be  developed. 

The  present  formulation  gives  separate  differential,  difference  and  non-linear  equations  for  the 
overall  network  and  for  each  of  the  individual  multiports.  This  leads  to  an  efficient  and  numerically 
stable  algorithm.  No  difficulties  have  been  encountered  in  analysing  networks  with  very  unequal 
time  constants  such  as  microwave  mixers. 

The  advantages  of  this  approach  are  summarized  below: 

(1)  A  large  network  can  be  divided  into  smaller  subnetworks  and  the  equations  for  each  subnet¬ 
work  are  derived  separately. 

(2)  A  library  of  subnetworks  can  be  developed  and  stored  for  future  use  without  the  need  of  an 
equivalent  circuit.  This  includes  transistors,  FETs,  diodes,  matching  sections,  filters  and 
couplers. 

(3)  The  equations  characterizing  a  non-linear  device  can  be  derived  to  match  experimental  data 
without  the  need  to  develop  a  physically  realizable  equivalent  circuit.  This  gives  a  greater 
flexibility  in  modelling  active  devices. 

(4)  The  subnetworks  developed  can  be  used  in  either  a  direct  integration  subroutine  or  a  harmonic 
balance  subroutine. 


2.  THE  GENERALIZED  L/D  MULTIPORT 

A  general  multiport  composed  of  individual  multiports  is  shown  in  Figure  1.  The  individual 
multiports  are  composed  of  lumped,  distributed  elements  and  dependent  sources.  The  lumped 
elements  are  linear  and  non-linear  resistors,  capacitors  and  inductors.  The  distributed  elements 
are  transmission  lines  coupled  or  uncoupled  embedded  in  homogeneous  or  inhomogeneous  media. 
The  overall  network  is  composed  of  all  individual  multiports  and  independent  sources.  Each 
multiport  has  current-driven  and  voltage-driven  ports.  These  are  ports  for  which  either  the  current 
or  the  voltage  is  considered  as  the  input.  The  /th  multiport  is  described  by 


x<  -  A>x>  +  B’ui  + 

(la) 

y>  =  C'x>  -F  D'u'  +  D'„u’„ 

(lb) 

Fi.  =  C[x'  +  D’u>  -F  D\„u'„ 

(Ic) 

where 

x>  =  (x{(/)x2(/)]^,x{{r)andx2(r)  are  the  lumped  and  distributed  state  vectors  of  the  /th 
multiport,  respectively. 


lodepeodent 

Qirrrat 

Sources 


lodependeDt 

Voltage 

Sources 


Figure  1.  General  lumped-distributed  non-linear  multiport 
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dx'(t)  „  ^  JT 

— ^x4(t-ET*) 


r*  is  the  delay  of  .he  ith  transmission  line. 


dt 

's  the  input  vector. 

WpiipP  is  the  output  vector. 

is  the  vector  of  the  controlling  voltages  and  currents  of  the  non-linear  elements, 
is  the  vector  of  the  non-linear  functions, 
are  real  matrices  of  the  state  and  output  equations  of  the  yth  multiport. 


ui 

y 

FL 

and  the  subscripts  cp  and  vp  refer  to  the  current-driven  or  voltage-driven  ports. 


With  each  non-linear  lumped  distributed  multiport  being  represented  by  equation  { 1 )  we  now 
proceeded  to  derive  the  state  equations  of  the  overall  network  which  consists  of  any  number  of 
individual  multiports. 


3.  THE  NETWORK  TOPOLOGY 

The  whole  network  is  obtained  by  interconnecting  all  multiports  and  independent  sources.  The 
topology  of  these  interconnections  is  represented  by  unconnected  graphs.  The  edges  of  each  graph 
have  to  satisfy  Kirchhoffs  laws.  A  forest  is  defined  and  Kirchhoffs  laws  can  be  expressed  in  the 
hybrid  form. 


it 

0 

D  ■ 

Vf 

.Vc  . 

_-D^ 

0  . 

Jc  . 

where  D  is  the  dynamical  transformation  matrix.'^ 


^vs.cp  1 

Ovs.vp 

1  ^\sxs 

D  = 

(2cp.cp : 

i 

^tp.vp  : 

^cp  .o» 

^vp.cp  1 

^vp.vp 

^vp.cs 

[^vs  icp. 

.1  ^Vp.f  j 

V,  = 

[Vv>V,p,, 

('^cp.cVvp.c'^rjl  ic  [icp.civp.cic^l 


(2) 


(3) 


and  the  subscripts  f,  c,  vs,  cs  refer  to  the  forest,  coforest,  independent  voltage  source  and 
independent  current  source  respectively. 

Let  us  define  the  following  vectors 


[I'cp.fVvp.c)  »  t/T  [Vvp.ft'cp.c]  ,  Ul, 

yi  [^cp.ft'vp.cl  ,  and  y2  [t'vp.fVcp.tl 

where  u,  is  the  source  vector  containing  all  the  independent  voltage  and  current  sources  of  the 
whole  network. 

It  should  be  noted  that  the  independent  voltage  and  current  source  edges  must  be  always  in 
the  forest  and  coforest,  respectively.  Without  loss  of  generality,  the  maximum  number  of  current- 
driven  ports  of  all  multiports  are  assigned  to  the  forest  (£>vp,cp  =  0). 

The  following  equations  can  be  obtained  from  (2)  and  (3). 


Ml  -  Fiyi  +  FiUj  +  Fim, 
>2  =  -Fiy,  +  FtU, 


{4a) 

(4b) 
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where 


F.  = 

0  ; 

^cp.vp 

.F,  = 

0 

1 

.-DJp.J 

0 

- 

0  . 

II 

r  0  ! 

I 

^cp.cs 

and  F4 

0  :  D..p,„  ■ 

_/)T  1 

L  ^vs.vp  ' 

0 

0  . 

Equations  (4)  are  auxiliary  equations  which  will  be  used  in  the  derivation  of  the  network  equations 
in  the  next  section. 


4.  FORMULATION  OF  THE  NETWORK  EQUATIONS 

The  state,  output  and  non-linear  equations  of  the  whole  network  consisting  of  a  number  of 
multiports  is  written  in  the  form. 


-f/»  -F  "t"  (^3) 

y,,  =  C,,x,,  +  D,,Ui,  +  D,„,u„  (5b) 

^nfy  {5c) 

where  ip,  Up,  u„  and  F„,,  are  real  vectors,  each  vector  contains  the  elements  of  the  corresponding 
vectors  of  all  multiports  (e.g.  Xp  =  in  is  the  number  of  all  multiports).  Ap.  Bp,  B„p, 

Cp.  Dp.  t)„p,  C’lp,  6,p,  and  D,„p  are  real  quasidiagonal  matrices,  each  matrix  contains  the  elements 
of  the  corresponding  matrices  of  all  multiports. 

TTie  state  vector  is  rearranged  to  contain  all  the  lumped  state  variables  followed  by  the  distributed 
ones  of  the  whole  network.  The  vectors  «p  and  y,,  arc  also  rearranged  according  to  the  forest  and 
coforest  edges  of  the  defined  vectors  u,.  ih.  >1  and  ys.  Hence  the  following  relations  are  obtained; 


.V  =  P,Xp 

(oa) 

u-  P.Up 

(6b) 

y  =  P^y,, 

(6c) 

where  F,  and  P2  are  elementary  transformation  matrices,  with  element  values  of  zero  or  one. 


u^[u^U2V,  y  =  (y,y2f 

From  (5)  and  (6)  the  overall  network  of  multiports  is  described  by 

x  =  ApX  +  BpU  +  B„pU„  (7a) 

y  =  C„.t  -t-  DpU  +  D„pU„  (7b) 

and 

F„  ==  CipX  +  0,pU  +  D,„pU„  (7c) 

where 


B„p  P 1 6„p 


Ap  =  P.ApPl, 


Bp  =  P,BpPi. 
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C,  =  PXpPJ.  Dp  =  D,p  =  P,[)„^ 

C.p  =  C,^PJ,  0,p  =  and  D,„p  =  D,„p 

Equation  (7b)  is  partitioned  as  follows. 


From  (4)  and  (8),  we  get 

iViiM  =  vi',jr  +  w,u,  +  IV, «„  (9) 


Table  I.  Stale-space  represenialton  ol  basic  lumped  elements 
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Table  H.  State-space  representation  of  transmission  lines 


Element 

Inputs 

Matn 

A 

ces  of  state  a 

B 

_ 

nd  output  equat 

C 

_ 

ions 

i 

D  ! 

; 

_ i 

— 

■ 

it' 

1-2 

0  -1 

-1  0 

■ 

1 

1 

■i 

n,  0 

0  Vo 

i 

u  = 

a 

.i2. 

1 

1 

0  Z„ 
Zo  0 

1 

B 

1 

Z„  0 

0  Zo 

il  i'2 

Transmisssion  line 

_ 

■ 

B 

1 

0  -r 
1-  oj 

1 

1 

0  1 
z.,  0 

1 

1 

0B2 

1 

1 

z„  0 

0  Yo 

! 

i 

u  ~ 

1 

1 

(1  1 

-1  1) 

1 

'o  z„ 

1  0 

1 

1 

-2K,  0 
0  2 

1 

1 

Y„  0 

0  Z„, 

,  Z„=l.'V',. 

-  -ilq - ^ 

11  =  fl 

0  1 

- 1  (1 

Cl 

1-:k,  ((| 

|T„! 

ii 

open-circuited  stub 

■ 

1 

(1 

z„ 

,,  Z„=1.K. 

-  ^ 

M  =  »'! 

0  -1 

-1  0 

0 

1, 

[-2V,,  01 

_ 1 

Short-circuited  stub 

u  -  il 

— 

0  -1 

1  0 

1) 

z„ 

i;  "1 

• 

(2n| 

where 


F 1  Dp  1 

-F,D, 

+  FI  Dp, 

FiDpz 

w,  = 


L-f/Cp,  -  c,.,  J 

-1 

F, 


F.j 


w,  = 


F,Dr. 


~FjD„pi  D„pzJ 


and  I„  is  a  unit  matrix. 

Finally,  the  network  equations  are  obtained  from  (7)  and  (9). 

X  -  Ax  +  Bu,  +  B„u„ 
y  =  Cx  +  Du,  +  D„u„ 


(10a) 

(10b) 
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Table  III.  State-space  representation  of  linear  controlled  sources 


Element 

Linear 

Excitation 

Matrices  of  state  and  output  equations 

relations 

inputs 

A 

B 

C 

D 

/T^ 

/2 

iH 

il  =  ail 

[il 

[0  0 

il 

in 

vl  =  0 

M  = 

v2 

a  0 

cccs 

il 

il 

- 1 

■ 

/jT', 

il  =  gil 

1 

vl 

[0  0 

vl 

T_- 

(1=0 

1 

v2 

g  oj 

vccs 

■ 

;1 

il 

A  * 

A 

1-2  =  (3vl 

'2] 

0  P' 

it 

?- 

C  =  0 

u  - 

.‘■J 

0  0. 

vevs 

il 

il 

4. 

n 

A 

v2  =  ril 

[ill 

bBi 

vl 

<^1 

vl  =  0 

M  = 

il 

1 

cevs 

F„  =  C,x  +  D,u,  +  Di^u„ 


(10c) 


A  =  flptvo'tv,  +  Ap 
B  =  BpWo'w2 
B„  =  flpWo 'iv,  + 

C=  D^iv,7'h',  +  Cp 
D  =  DpWii  ^W2 
D„  =  DpW^'wy  +  D„p 

=  DipW„'»v,  +  C|p 
Oi  =  DipWo ‘w^ 

Dx„  =  DipW^ ‘wi  -1-  D,„p 


where 
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Table  IV.  State-space  representation  of  basic  non-iinear  elements 


Element 


Voltage-controlled 

resistor 


Current-controlled 

resistor 


C\ 


1 


Non-linear  capacitor 


j  -=:|} 


Non-linear  inductor 


Non-linear  voltage 
source 


Non-linear  current 
source 


Non-linear 

relations 


f  =  /(V) 


V  =  f(i) 


>•  =  g(<?) 

E  =  g(v„.  C.)-v<, 


<  =  /(4<) 
7  —  /(fo, 


E  = 


J  =  f{x,u.u„.i) 
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Figure  2.  Circuit  diagram  and  multiport  equivalent  of  balanced  mixer,  (a)  Schematic  diagram  of  balanced  mixer,  (b) 

Multiport  equivalent  of  balanced  mixer 


The  matrix  w„  may  be  singular  due  to  the  dependence  between  some  of  the  lumped  state  variables. 
Such  dependence  which  is  due  to  the  interconnection  of  all  multiports  can  only  arise  under  the 
following  conditions; 

( 1 )  The  network  has  some  cutsets  consisting  of  only  inductors  and  current  sources. 

(2)  The  network  has  some  loops  consisting  of  only  capacitors  and  voltage  sources. 

(3)  The  presence  of  dependent  sources  in  some  special  cases.  This  condition  does  not  occur  in 
practical  networks. 

The  dependent  state  variables  can  be  eliminated  by  elementary  row  and  column  operations  on 
the  coefficient  matrices  in  (9). 


5.  SIMULATION 

A  general  computer  program  has  been  developed  for  the  analysis  of  non-linear  L/D  networks. 
The  formulation  of  the  network  equations  has  been  established  by  using  sparse  matrix  techniques. 
The  solution  of  (10)  can  be  obtained  as  explained  in  Reference  1. 

The  explicit  forms  of  the  matrices  of  network  equations,  describing  the  devices  and  subcircuits 
commonly  used  are  implemented  in  the  program.  The  advantage  of  the  proposed  method  is  that 
the  developed  program  deals  with  these  circuits  as  multiports,  describing  their  terminal  behaviour 
instead  of  dealing  with  their  basic  circuit  elements.  Basic  linear  and  non-linear  circuit  elements 
(such  as  resistors,  inductors,  capacitors,  controlled  sources  and  transmission  lines)  can  also  be 
represented  as  multiports.  The  state  space  representation  of  these  elements  is  given  in  Tables  I. 
II,  III  and  IV, 


6.  EXAMPLES 

The  developed  program  has  been  applied  to  several  examples.  In  the  following  examples,  the 
circuit  is  partitioned  into  multiports  using  some  of  the  implemented  subcircuits  in  the  program 
such  as  diodes,  MESFET’s  and  matching  sections. 
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6. 1.  Schottky  diode  balanced  mixer 

A  balanced  microwave  mixer  circuit  using  two  silicon  Schottky  diodes  DC1533G  was  analysed. 
The  local  oscillator  and  intermediate  frequencies  are  7  and  0144  GHz.  respectively.  The  schematic 
diagram  of  the  mixer  is  shown  in  Figure  2(a).  The  network  is  divided  into  a  number  of  multipons 
and  models  of  each  multiport,  including  the  Schottky  diodes,  are  developed  and  scoicd  in  the 
program  library.  The  overall  network  is  then  analysed  as  an  interconnection  of  the  muliiports.  as 
shown  in  Figure  2(b).  The  equivalent  representations  of  each  multiport  are  given  in  Table  V.  A 
higher  level  of  hierarchy  is  also  possible  and  larger  multiport  representations  can  be  made  if 
required.  The  output  waveforms  before  and  after  IF  filter  are  snown  in  Figure  3.  The  variation 
of  the  conversion  loss  with  RF  frequency  is  shown  in  Figure  4. 


6.2.  MESFET  frequency  doubler 

A  similar  procedure  has  been  used  to  analyse  a  2-5  GHz  frequency  doubler,  using  a  Plessey 
P35-1105-1  MESFET,  shown  in  Figure  5(a)  and  the  multiport  equivalent  is  shown  m  Figure  5(b). 
The  output  waveform  is  shown  in  Figure  6.  The  output  is  further  analysed  and  the  fr;quency 
response  is  obtained.  The  circuit  has  been  built  and  tested  and  the  theoretical  frequency  response 
is  compared  with  the  practical  measurements  in  Figure  7.  Good  agreement  is  shown  between 
measured  and  predicted  results  which  gives  confidence  in  the  developed  method. 


7.  CONCLUSION 

Non-linear  lumped-distributed  networks  can  now  be  analysed  in  the  time  domain  as  an  intercon¬ 
nection  of  multiports.  The  overall  network  is  divided  into  a  number  of  subnetworks  and  each 


Figure  3.  Simulated  output  of  microwave  mixer 


initscutt  8 
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Figure  4.  Frequency  response  of  microwave  mixer 


T 


(a)  Schematic  diagram  of  frequency  doubler 


X  (b)  Multiport  equivalent  of  frequency  doubler 


Figure  5.  Circuit  diagram  and  multiport  equivalent  of  frequency  doubler,  (a)  Schemattc  diagram  of  frequency  doubler, 
(b)  Circuit  diagram  and  multiport  equivalent  of  frequency  doubler 


subnetwork  is  characterized  separately.  A  library  of  subnetworks  can  be  developed  from  active 
elements  such  as  transistors,  FETs,  diodes,  eU.  with  very  little  storage  required.  The  non-linear 
multiports  can  be  used  in  either  a  direct  integration  subroutine  or  using  the  harmonic  balance 
method. 
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Figure  (>.  Output  waveform  of  frequency  doubler 


tccond  HanMOic 


Figure  7.  Frequency  response  of  frequency  doubler 
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